the weather impact on robusta coffee yield in Vietnam

Keywords: modélisation statistique, régression et classification, estimation du rendements agricoles, validation croisée, calibration du modèle climatique, impact du changement climatique statistical modelling, regression and classification, crop yield estimation, cross-validation, calibration of climate model, climate change impact. vii Modélisation statistique, régression et classification, estimation du rendements agricoles, validation croisée, calibration du modèle climatique, impact du changement climatique Statistical modelling, regression and classification, crop yield estimation, cross-validation, calibration of climate model, climate change impact

Weather and climate strongly impact crop yields. Many studies based on different techniques have been done to measure this impact. This thesis focuses on statistical models to measure the sensitivity of crops to weather conditions based on historical records. When using a statistical model, a critical difficulty arises when data is scarce, which is often the case with statistical crop modelling. There is a high risk of overfitting if the model development is not done carefully. Thus, careful validation and selection of statistical models are major concerns of this thesis.

Two statistical approaches are developed. The first one uses linear regression with regularization and leave-one-out cross-validation (or LOO), applied to Robusta coffee in the main coffee-producing area of Vietnam (i.e. the Central Highlands).

Coffee is a valuable commodity crop, sensitive to weather, and has a very complex phenology due to its perennial nature. Results suggest that precipitation and temperature information can be used to forecast the yield anomaly with 3-6 months' anticipation depending on the location. Estimates of Robusta yield at the end of the season show that weather explains up to 36 % of historical yield anomalies.

The first approach using LOO is widely used in the literature; however, it can be misused for many reasons: it is technical, misinterpreted, and requires experience.

As an alternative, the "leave-two-out nested cross-validation" (or LTO) approach, is proposed to choose the suitable model and assess its true generalization ability.

This method is sophisticated but straightforward; its benefits are demonstrated for Robusta coffee in Vietnam and grain maize in France. In both cases, a simpler model with fewer potential predictors and inputs is more appropriate. Using only the LOO method, without any regularization, can be highly misleading as it encourages choosing a model that overfits the data in an indirect way. The LTO approach is also useful in seasonal forecasting applications. The end-of-season grain maize yield estimates suggest that weather can account for more than 40 % of the variability in yield anomaly. Climate change's impacts on coffee production in Brazil and Vietnam are also vi studied using climate simulations and suitability models. Climate data are, however, biased compared to the real-world climate. Therefore, many "bias correction" methods (called here instead "calibration") have been introduced to correct these biases. An up-to-date review of the available methods is provided to better understand each method's assumptions, properties, and applicative purposes. The climate simulations are then calibrated by a quantile-based method before being used in the suitability models. The suitability models are developed based on census data of coffee areas, and potential climate variables are based on a review of previous studies using impact models for coffee and expert recommendations. Results show that suitable Arabica areas in Brazil could decrease by about 26 % by the mid-century in the high-emissions scenario, while the decrease is surprisingly high for Vietnamese Robusta coffee (≈ 60 %). Impacts are significant at low elevations for both coffee types, suggesting potential shifts in production to higher locations.

The used statistical approaches, especially the LTO technique, can contribute to the development of crop modelling. They can be applied to a complex perennial crop like coffee or more industrialized annual crops like grain maize. They can be used in seasonal forecasts or end-of-season estimations, which are helpful in crop management and monitoring. Estimating the future crop suitability helps to anticipate the consequences of climate change on the agricultural system and to define adaptation or mitigation strategies. Methodologies used in this thesis can be easily generalized to other cultures and regions worldwide.

Résumé

La météo et le climat ont un impact important sur les rendements agricoles. De nombreuses études basées sur différentes approches ont été réalisées pour mesurer cet impact. Cette thèse se concentre sur l'utilisation de modèles statistiques pour mesurer la sensibilité des cultures aux conditions météorologiques sur la base des enregistrements historiques. Lors du développement et de l'utilisation d'un modèle statistique, une difficulté critique survient lorsque les données sont rares, ce qui est souvent le cas pour la modélisation des cultures. Il y a un risque élevé de sur-apprentissage si le modèle n'est pas développé avec certaine précautions.

Ainsi, la validation et le choix du modèle sont deux préoccupations majeures de cette thèse.

Deux approches statistiques sont développées. La première utilise la régression linéaire avec régularisation et validation croisée (c.-à.-d. leave-one-out ou LOO), appliquée au café robusta dans la principale région productrice de café du Vietnam (c.-à.-d. les Montagnes Centrales). Le café est une culture rémunératrice, sensible aux intempéries, et qui a une phénologie très complexe en raison de sa nature pérenne. Les résultats suggèrent que les informations sur les précipitations et la température peuvent être utilisées pour prévoir l'anomalie de rendement avec une anticipation de 3 à 6 mois selon la région. Les estimations du rendement du robusta à la fin de la saison montrent que les conditions météorologiques expliquent jusqu'à 36 % des anomalies de rendement historiques. Cette première approche de validation par LOO est largement utilisée dans la littérature ; cependant, elle peut être mal utilisé pour de nombreuses raisons : elle est technique, mal interprétée et nécessite de l'expérience. Une alternative, l'approche "leavetwo-out nested cross-validation" (ou LTO), est proposée pour choisir le modèle approprié, évaluer sa véritable capacité de généralisation et choisir la complexité du modèle optimale. Cette méthode est sophistiquée mais simple. Nous démontrons son applicabilité pour le café robusta au Vietnam et le maïs en France.

Dans les deux cas, un modèle plus simple avec moins de prédicteurs potentiels et d'entrées est plus approprié. Utiliser uniquement la méthode LOO, sans aucune iv régularisation, peut être très trompeur car cela encourage à choisir un modèle qui sur-apprend les données de manière indirecte. L'approche LTO est également utile dans les applications de prévision saisonnière. Les estimations de rendement du maïs en fin de saison suggèrent que les conditions météorologiques peuvent expliquer plus de 40 % de la variabilité de l'anomalies de rendement en France.

Les impacts du changement climatique sur la production de café au Brésil et au Vietnam sont également étudiés à l'aide de simulations climatiques et de modèles d'impact (appelées "suitability models"). Les données climatiques sont cependant biaisées par rapport au climat réel. Par conséquent, de nombreuses méthodes de "correction de biais" (appelées ici "calibration") ont été introduites pour corriger ces biais. Une présentation critique et detaillée de ces calibrations dans la littérature est fournie pour mieux comprendre les hypothèses, les propriétés et

Acknowledgment

I did not realize three years went by so quickly until I sat down to write these lines. The last three years have been valuable for me to learn and grow in many aspects, professionally and personally. This work would not have been possible without the support of many people: my supervisor, thesis committee, colleagues, administrative team, friends, and family. I thank all for their companionship, assistance, and motivation. I also thank the readers for their interest in this research.

Thầy Filipe, saying thank you is never enough to capture my gratitude for everything you have done for me. I am so lucky to have you as my supervisor who always guides and encourages me. You have taught me so much, both in science and in life! It is my great pleasure to have Prof. Dr. Christoph Gornott and Assoc. Prof. Erwan Personne as reporters on my jury. I also thank Dr. Frédérique Cheruy, Dr. Tamara Ben Ari, Prof. Dr. Philippe Ciais, and Dr. Eric Rahn, for accepting my invitation to become a part of the jury.

Many thanks to Eric Rahn and Pablo Imbach for opening the door to CIAT and the coffee world and Philippe Ciais for opening the door to CLAND. Eric, thanks so much for your continued support, comments, and collaboration with the coffee work. My thesis could not be finish without your help. I acknowledge the fruitful scientific discussions during the two thesis monitoring committees with Dr. David Makowski, Dr. Benoit Bertrand, and Dr. Nathalie de Noblet. I am very grateful for their comments, suggestions, and advice. I also highly appreciate the comments and suggestions from Eric Rahn, Victor Pellet, and Alain Coulais. Their comments have been very thorough and helpful in improving the thesis manuscript. I would like to thank Tarikakan Demirturk x for his IT support during my three years at LERMA. Special thanks to Catherine Prigent and Alain Coulais for their kindness, friendship, and encouragement. Also, many thanks to all my colleagues and friends at Paris Observatory, especially Nguyen Thu Hang, Eulalie Boucher, Matthew Heberger, Victor Pellet, Qinwei Ran, Elise Blanchard, Murielle Chevrier and others for all the time we have shared in the last three years.

Most importantly, I owe my family and friends who have always been there for me. My parents always support me unconditionally and give me the freedom to discover the world. Special thanks to my beloved friend, Đỗ Huy Duy, who has always been with me, supported and cared for me.

Finally, I thank myself for not giving up on me. There were hard times and even depression (for many personal reasons), but "we" made it to the end. I dedicate it to "you".

xi "Climate change impacts are stressing agriculture, forestry, fisheries and aquaculture, increasingly hindering efforts to meet human needs." [START_REF] Kerr | Food, fibre, and other ecosystem products[END_REF] 1.1 General background and context

What is agriculture and why is it important?

Agriculture is the basic material production of society, including all forms of activities connected with growing, harvesting and preliminarily processing crops; breeding, raising and caring for animals; and tending gardens and nurseries [ILO, 1999]. The history of agriculture began thousands of years ago. For instance, wild grains were collected and eaten at least 105 000 years ago [START_REF] Harmon | Humans feasting on grains for at least 100,000 years[END_REF]. Cereal cultivation of emmer, barley, and oats was observed near the sea of Galilee about 23 000 years ago [START_REF] Snir | The origin of cultivation and proto-weeds, long before neolithic farming[END_REF]. Animals such as cattle, sheep, and pigs 1

Chapter 1 Introduction were domesticated over 10 000 years ago [START_REF] Larson | Worldwide phylogeography of wild boar reveals multiple centers of pig domestication[END_REF][START_REF] Mctavish | New world cattle show ancestry from multiple independent domestication events[END_REF]. For decades, agriculture involved the production of important food crops.

Although humanity has changed significantly, agriculture remains and has become increasingly important. Many raw materials (such as rubber, cotton, sugar, or oil) come from agriculture. These materials are essential to major industries in many ways (e.g. the manufacturing of pharmaceuticals, diesel fuel, fibres, and more). In addition, raw materials from agriculture constitute a large part of what is traded internationally. In 2021, the European Union (EU) maintains its position as the top trader in the world's agricultural products, both in terms of exports and imports [Directorate-General for Agriculture and Rural Development, 2022].

Agriculture is also an important economic sector in the economy of many countries. For instance, according to the Food and Agriculture Organization (FAO), agriculture accounted for 4 % of global gross domestic product (GDP). In some developing countries, it accounted for more than 25 % of GDP in 2018. Yet, with the challenges in today's world -e.g. population growth, Covid-19 pandemic, conflicts, crises, and natural disasters -we need to draw tremendous attention to the importance of this sector. This thesis thus focuses on one of the fundamental branches of agriculture -crop yield/production.

How is crop yield variability explained?

In agriculture, crop yield is strongly influenced by several factors [START_REF] Liliane | Factors affecting yield of crops[END_REF] such as technology (e.g. agricultural practices, managerial decisions), biology (e.g. genetics, diseases, pests, weeds) and the environment (e.g. climatic conditions, soil fertility, topography, water quality). These factors can pose a significant risk to farms leading to yield reduction when they are not correctly monitored and well managed. They often interact with each other and directly impact the crop yield variability; some factors have a short-term effect while others act over the long term. Understanding different sources of yield variability can support farmers in adapting their agricultural system to be more resilient. Generally, the variability of meteorological conditions within the limits of usual climatic conditions explains most of the annual crop yield variability over a short period. On the other hand, agricultural practices and/or the introduction of new varieties (i.e. biological factors) account for most of the variability over a long-term period, often in the form of a trend over ten to twenty or thirty years [START_REF] Schauberger | Yield trends, variability and stagnation analysis of major crops in France over more than a century[END_REF]. For extended periods, climate change and environmental improvement or degradation (e.g. in soil properties and water quality) might well be other important factors influencing yield [START_REF] Bouman | Yield prediction by crop modelling/remote sensing[END_REF][START_REF] Mueller | Closing yield gaps through nutrient and water management[END_REF]. These various elements are not independent of each other. For example, agricultural practices can be partly dependent on weather conditions: certain practices diminish or, in contrast, reinforce the variability due to weather conditions. Therefore, the effects of weather and those of agricultural practices (and/or other factors) cannot always be easily distinguished. In this study, the long-term evolution will be considered in the form of a trend, making it possible to introduce the inter-annual variability of crop yields: the weather remains the major uncontrollable factor influencing the development of crops [START_REF] Taylor | Weather and yield trends[END_REF][START_REF] Ray | Climate variation explains a third of global crop yield variability[END_REF]. The weather can vary considerably from field to field and from year to year. Thus, understanding the relationship between weather variability and crop production can help assess agricultural production systems' resilience to future climate conditions and help identify appropriate adaptation strategies to climate change.

Weather and climate impact on crop yield

Many studies have shown the importance of weather and climate in explaining part of crop yields. Weather can affect crop production in all phases of the crop growing cycle. For example, unfavourable weather conditions in the early growth stages of grain maize may limit the size of the leaves and thus the photosynthetic capacity [START_REF] Darby | Plant physiology: Critical stages in the life of a corn plant[END_REF]. In later stages, adverse conditions (heatwave and drought) may reduce the number of silks produced, resulting in poor pollination of the ovules and restricting the number of kernels that develop; or growth may stop prematurely and restrict the size of the kernels produced [START_REF] Ritchie | How a corn plant develops[END_REF]. Another example is a study that analysed weather impacts on Robusta coffee yield in Vietnam and Indonesia suggesting a 14 %yield reduction for every 1 • C increase in mean minimum/maximum temperatures (above 16.2/24.1 • C) during the growing season [START_REF] Kath | Not so robust: Robusta coffee production is highly sensitive to temperature[END_REF]. Notably, extreme weather events (frosts, heat waves, or prolonged drought) have severe detrimental effects on plant growth and development, and thus, crop production [START_REF] Hatfield | Temperature extremes: Effect on plant growth and development[END_REF][START_REF] Mathieu | Using Neural Network Classifier Approach for Statistically Forecasting Extreme Corn Yield Losses in Eastern United States[END_REF][START_REF] Vogel | Identifying meteorological drivers of extreme impacts: An application to simulated crop yields[END_REF]. On a larger scale, the study of [START_REF] Ray | Climate variation explains a third of global crop yield variability[END_REF] suggested that climate variability accounts for roughly a third of the observed yield variability of maize, rice, wheat and soybean worldwide. This value can even be >60 % for maize in some specific regions, e.g. in numerous counties of the U.S. and some provinces in China. Due to the importance of weather-to-yield relation, research on crop yield simulation -which can cover several tasks such as end-of-season yield estimations, seasonal yield forecasts, or long-term predictions -is essential. Agricultural yield estimations (i.e. the real-time quantification of yield potential) are useful for the whole agricultural production chain, from the farmers, traders or insurers to policymakers. For instance, farmers can adapt their management, and traders or insurers can start planning their pricing. Governments can prepare their food balance reports used for market analyses and decisions related to managing stocks, imports, and exports. Reliable yield simulation (e.g. simulation/forecasts within the current growing season) can be used strategically to analyse the possible consequences of management practices. It can also be used as a tactical management tool, e.g. sowing or harvesting date. If the predictors are made in real-time, they can be used directly by farmers to plan production (i.e. quick decision-making).

When it comes to climate change, agriculture is among the most vulnerable sectors. The impacts of climate change on agriculture are unevenly distributed around the world, depending on crops, regions, timeframe, and projection scenarios. Although positive effects have been identified in several regions (e.g. in the high latitudes [START_REF] Porter | Food security and food production systems[END_REF]), climate change will disproportionately affect crop yields among regions, with more negative than positive effects being expected in most areas [START_REF] Pörtner | Climate Change 2022: Impacts, Adaptation, and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF]. According to the Intergovernmental Panel on Climate Change (IPCC) sixth assessment report [START_REF] Kerr | Food, fibre, and other ecosystem products[END_REF], the overall median per decade effect is -2.3 % for maize, -3.3 % for soybean, -0.7 % for rice and -1.3 % for wheat. [START_REF] Chen | Impacts of climate change on agriculture: Evidence from china[END_REF] stated that global warming had caused an economic loss of about $820 million to corn and soybean producers in China in the past 10 years. In another example, Bunn et al. [2015b] showed that climate change could reduce the global area suitable for coffee by about 50 % across emissions scenarios. Furthermore, climate change is likely to increase the risk of food insecurity, especially in developing countries [FAO, 2015]. Adaptation options for crop production systems are assessed using crop models together with climate simulations. These models can then be used to project possible future crop yields with respect to the change in climate conditions. Such projections can support farmers' decision-making in their management practices (e.g. changing the use of fertilisers or even acting on irrigation), the modification of crop systems (e.g. changing crop varieties, adjusting sowing dates), the extent of crop propagation, and/or other adaptation strategies.

Methodological approaches used in crop modelling

Several modelling approaches have been used to study the relationship between weather/climate and crop yield. The two main approaches are processbased models 1 [van Oijen et al., 2010; Dormann et al., 2012; Rahn et al., 2018; Vezy et al., 2020] and statistical models [START_REF] Lobell | On the use of statistical models to predict crop yield responses to climate change[END_REF][START_REF] Ray | Climate variation explains a third of global crop yield variability[END_REF][START_REF] Gornott | Statistical regression models for assessing climate impacts on crop yields: A validation study for winter wheat and silage maize in Germany[END_REF], even though hybrid approaches have also been developed (e.g. functional model in [START_REF] Chipanshi | Evaluation of the integrated canadian crop yield forecaster (iccyf) model for in-season prediction of crop yield across the canadian agricultural landscape[END_REF]). Although the approaches overlap on some points, each has its own advantages and disadvantages, especially in estimating crop yields.

Process-based models use mathematical equations to simulate the mechanisms that control crop growth and/or seed formation based on soil characteristics, weather, initial conditions, and management practices [START_REF] Hoogenboom | From genome to crop: integration through simulation modeling[END_REF][START_REF] Van Oijen | Coffee agroforestry systems in Central America: II. Development of a simple processbased model and preliminary results[END_REF]. This type of model is also a powerful tool for crop yield predictions, particularly at a field scale, as it provides a clear physiological mechanism for linking weather to crop yield outcome [START_REF] Roberts | Comparing and combining process-based crop models and statistical models with some implications for climate change[END_REF]. Once properly calibrated and evaluated with observed data, a process-based model can also provide an understanding of more profound physiological constraints and thus help the targeted development of adaptive technologies; or it can provide information on possible management interventions to better cope with climate change [START_REF] Laux | Impact of climate change on agricultural productivity under rainfed conditions in cameroon-a method to improve attainable crop yields by planting date adaptations[END_REF]. However, as this model is originally developed for local field conditions, large-scale implementations might need an extensive set of input data, such as soil information, meteorological variables, eco-physiological parameters describing the crop variety, and agro-management practices (planting and harvesting dates, irrigation, and fertilisation). Obtaining these datasets at a large scale (i.e. regional or global scales) is not easy. Another challenge with process models is that they are typically parameterized using data from experimental field plots; therefore, their use in other regions would require a re-calibration [START_REF] Folberth | Regionalization of a large-scale crop growth model for sub-saharan africa: Model setup, evaluation, and estimation of maize yields[END_REF]. To this end, due to intensive data and calibration requirements of process-based crop models, the use of detailed process-based crop models for timely predictions of crop yield at regional or global scales remains challenging [START_REF] Lobell | On the use of statistical models to predict crop yield responses to climate change[END_REF].

As an alternative to process-based models, statistical models tend to operate at larger scales, often compatible with the production statistics at administrative levels collected by national governments. This approach has been widely used over the last decade (see, for example, [START_REF] Lobell | On the use of statistical models to predict crop yield responses to climate change[END_REF]; [START_REF] Gornott | Statistical regression models for assessing climate impacts on crop yields: A validation study for winter wheat and silage maize in Germany[END_REF]; [START_REF] Mathieu | Statistical weather-impact models: An application of neural networks and mixed effects for corn production over the United States[END_REF]; [START_REF] Kern | Statistical modelling of crop yield in central europe using climate data and remote sensing vegetation indices[END_REF]), and the aim is to find relations between a set of explanatory variables and crop yield. The explanatory variables are either only weather information or a combination of weather and non-weather factors. In contrast to the process-based approach, this modelling process -also called the data-driven approach -starts with the data and aims to extract as much information from it as possible, with the fewest number of assumptions. Therefore, statistical models often have several advantages, such as easy calculations and less time-consuming. In addition, statistical approaches usually require fewer input data, thus, facilitating their application and scalability to other locations with much more simplicity than process-based models. However, by not including the plant-physiological interactions (e.g. biomass growth, photosynthesis, nutrient uptake), statistical models are thus limited in the information they can convey to the farmers, e.g. whether or not to add fertilisers. These models can also be misused, and they are limited in the information they can provide outside the range of considered input variables. Also, statistical models do not consider the complex and continuous soil-plant-atmosphere system, which can be critical when dealing with locations with very different soils. For instance, the response of crops to a certain amount of rainfall on loam soil is different from that on clay soil. Nevertheless, statistical models can take into account the impacts of various factors and capture their differences, e.g. water and/or heat stress during the flowering period might reduce yield more than the stress during the vegetative phase. Besides the weather impacts, statistical models can also consider the impacts of socio-economic factors (indirect factors), for example, in [START_REF] Gay | Potential impacts of climate change on agriculture: A case of study of coffee production in veracruz, mexico[END_REF]; [START_REF] Conradt | Extending and improving regionalized winter wheat and silage maize yield regression models for germany: Enhancing the predictive skill by panel definition through cluster analysis[END_REF]. Time-invariant yield impacts [START_REF] Hoffman | Analysis of climate signals in the crop yield record of sub-saharan africa[END_REF] can be included as well, for instance, in the constant term of linear regression (intercept).

Overall, by including many fundamental processes in the plant life cycle, a detailed process-based model can probably better explain the link between weather and crop yields. However, this model type is more appropriate at a local scale when much information is available to describe all physical processes. On the other hand, a statistical model aims only at explaining a part of the yield variability (e.g. the part that is driven by weather conditions). Consequently, less information is required for a less ambitious target. This statistical approach thus can leverage the available data to be scalable, parsimonious, and generalizable on a global scale.

The research focuses 1.2.1 Difficulties in statistical crop modelling

In this thesis, we first focus on developing statistical crop models to study the weather's impact on crop yields. A very significant difficulty for these models comes from the limited amount of data in historical databases (or from a large number of missing values). In agriculture, one year of (yield) data represents one sample, and it is often difficult to obtain data for more than two or three decades. For instance, [START_REF] Prasad | Crop yield estimation model for Iowa using remote sensing and surface parameters[END_REF] built crop yield estimation models for corn and soybean with 19 years of yield data. [START_REF] Ceglar | Impact of meteorological drivers on regional inter-annual crop yield variability in France[END_REF] studied the impact of meteorological drivers over 26 years on grain maize and winter wheat yield in France. A small sample size poses two challenges to crop modellers. First, it makes it hard to choose among competing models. Second, it is difficult to assess the quality of the chosen model. To overcome these small sample size issues, cross-validation is used in many crop models [START_REF] Kogan | Winter wheat yield forecasting in Ukraine based on Earth observation, meteorologicaldata and biophysical models[END_REF][START_REF] Zhao | Comparing empirical and survey-based yield forecasts in a dryland agro-ecosystem[END_REF][START_REF] Li | Toward building a transparent statistical model for improving crop yield prediction: Modeling rainfed corn in the u.s[END_REF]. A common type of cross-validation is the leave-one-out crossvalidation (LOO). However, this technique can sometimes become a bad practice, especially with a limited number of samples, by suggesting an overly complex model (see, for instance, in [START_REF] Jayakumar | Climate-based statistical regression models for crop yield forecasting of coffee in humid tropical Kerala, India[END_REF]; de Oliveira Aparecido et al. [2017]; [START_REF] Niedbała | Application of multiple linear regression for multi-criteria yield prediction of winter wheat[END_REF]). Careful selection and validation of statistical models are thus major concerns of this thesis. We aim to find the right statistical approaches and contribute to developing a new modelling methodology.

Calibrating climate change simulations

We will also perform climate change impact assessments with the example of coffee. To accomplish this task, climate models are often used as the essential inputs providing information about climate change. However, these climate models often suffer from substantial biases and errors compared to the current observations on which the impact models are calibrated. Therefore, many users of climate model data apply some form of "bias correction" (called here instead "calibration" techniques) to calibrate the climate simulations towards the observations. Several studies have attempted to review them, for example, the reviews of [START_REF] Teutschbein | Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods[END_REF]; [START_REF] Maraun | Bias correcting climate change simulations -a critical review[END_REF] or the comparisons in [START_REF] Watanabe | Intercomparison of bias-correction methods for monthly temperature and precipitation simulated by multiple climate models[END_REF]; [START_REF] Gutjahr | Comparing precipitation bias correction methods for high-resolution regional climate simulations using cosmo-clm[END_REF]. Since then, numerous novel methods have been proposed to correct and minimise climate model errors. For instance, some studies may aim at adjusting the wet-day frequencies, while others may prefer to consider the recurrence interval value. This variety is indeed an advantage as more choices allow the climate modellers to choose the best one for their own particular needs. Nonetheless, the variety of solutions can also become a difficulty when the impact modellers do not have the statistical expertise to make a choice among them. Consequently, an up-to-date and synthetic review of existing methods is necessary and will be one of the main focuses of this thesis.

Case studies

To further investigate and understand the impacts of weather and climate conditions on crop yields, we will focus on a very challenging crop, which is coffee. Unlike other common crops, coffee plants are woody perennial evergreens that can be productive for more than 50 years, and they often demonstrate a very complex two-year phenological cycle [Wintgens, 2004a;[START_REF] Valadares | Yield and production bienniality of dense coffee plantations under different levels of N and K[END_REF][START_REF] Aparecido | Forecasting of the annual yield of Arabic coffee using water deficiency[END_REF]. Furthermore, the management practices vary significantly among coffee-growing countries, even for the two leading producers (i.e. Brazil and Vietnam). In addition, coffee yield data, especially at a regional Chapter 1 Introduction scale, are often not freely accessible. Despite several challenging factors, we selected this crop for several reasons. First, coffee is one of the world's most valuable and widely traded commodity crops. It is a major cash crop and an important foreign exchange for many countries, including those in Latin America and the Caribbean, South and Southeast Asia, and Africa. Second, coffee has proven to be highly sensitive to weather [START_REF] Craparo | Coffea arabica yields decline in tanzania due to climate change: Global implications[END_REF][START_REF] Kath | Not so robust: Robusta coffee production is highly sensitive to temperature[END_REF][START_REF] Kath | Temperature and rainfall impacts on robusta coffee bean characteristics[END_REF] and climate [START_REF] Davis | The Impact of Climate Change on Indigenous Arabica Coffee (Coffea arabica): Predicting Future Trends and Identifying Priorities[END_REF]Bunn et al., 2015b;[START_REF] Moat | Resilience potential of the ethiopian coffee sector under climate change[END_REF]. Although several impact studies have addressed the interactive effects of weather and climate conditions on coffee plantations, few of them have access to the long-term coffee yield time series at regional scales. They often worked on some specific locations (e.g. few Brazilian municipalities or farm-level data) with limited time series (the maximum of 10 years) [de Oliveira Aparecido et al., 2017;[START_REF] Valeriano | Estimation of coffee yield from gridded weather data[END_REF][START_REF] Kath | Not so robust: Robusta coffee production is highly sensitive to temperature[END_REF]. Some other studies made use of indirect coffee information, such as the occurrence of coffee production without yield data [Bunn et al., 2015a,b]. Finally, several studies have discussed climate change adaptation of coffee production at local and even global scales [START_REF] Läderach | Climate change adaptation of coffee production in space and time[END_REF]Bunn et al., 2015b]. However, there is still room for improvement concerning: the coffee input database, the models used, and the selection and calibration of climate simulations. In addition, an application to French grain maize is also presented to show that the methodologies proposed here can easily be extended to other regions and crops (in particular, crops with simpler phenology).

Objectives and plan of the thesis

This thesis aims to improve the statistical modelling of crop yield, especially coffee. The statistical approach is chosen here because it is data-driven and, therefore, can be more easily extended globally. The developed techniques are then applied to coffee production to better understand weather and climate's impacts on agricultural yields. To achieve these aims, the thesis covers three main objectives.

• First, we measure the sensitivity of crops to weather conditions, focusing on two main types of coffee in the two leading countries, i.e. Arabica coffee in Brazil and Robusta coffee in Vietnam.

• Then, we focus on developing statistical crop yield models. Two approaches are performed: (1) a common approach using a simple regression model with regularization techniques and careful validation, and (2) a proposed approach using statistical models with a robust validation tool (which is used to choose the right model and assess its quality). Both approaches will be first applied to Robusta coffee as it is highly sensitive to weather. The

Chapter 1 Introduction first one is widely used in the literature; however, it is often technical, requires experience, and sometimes is challenging to interpret, especially for non-statisticians. The second -proposed -approach is more advanced but is easier to implement. This direct implementation will be shown not only in the Robusta coffee application but also in a more industrialized crop like French grain maize. Comparisons are made between our proposed approach and the common (but can be misleading) cross-validation approach.

• Finally, we investigate the climate change impact on coffee production. As for the long-term assessments, many factors can change (significantly), including weather conditions, management practices, or agricultural policies.

Estimating the future yield values might not be practical in this context. Therefore, we will instead present a suitability assessment, showing how climate change will impact the coffee production area. In addition, a synthetic review of available methods to calibrate climate model data is introduced.

After this introduction (chapter 1), the thesis is structured into seven chapters. Chapter 2 introduces the datasets, including agricultural and weather/climate data. Although we only focus on two crops, many available datasets are introduced for future studies. Chapter 3 addresses the methodology from a mathematical point of view (for statistical modelling) and a modelling point of view. Different regularization and especially cross-validation techniques are discussed here. Chapter 4 focuses on Robusta coffee yield estimation in Vietnam. Although several analyses are done for both Arabica coffee and Robusta coffee, Arabica coffee has ultimately been excluded due to the coarse resolution of the available yield data. Two statistical impact models are then developed and compared for Robusta coffee. This is followed by additional applications using the second approach on French grain maize, presented in chapter 5. Chapter 6 reviews existing calibration methods, and then chapter 7 describes the climate change assessments on coffee production in Brazil and Vietnam. Finally, chapter 8 concludes this thesis and provides scientific perspectives for future works. In this chapter, we will go through the available datasets that (potentially) will be used in the thesis. First, the agricultural data, including the coffee and European Union (EU) crops databases, are described. The coffee database introduces an overview of global coffee cultivation, the administrative division, and the available coffee databases worldwide. Here, we will focus more on the available data of the two leading coffee-producing countries (i.e. Brazil and Vietnam). Some other coffee-related data, which are needed for our impact assessments (in chapter 7), are also presented. In the case of the EU crops databases, there are interesting crop data at sub-national levels for many EU countries with high quality; however, this thesis focuses mainly on French crops (i.e. grain maize) due to time limitations. Nevertheless, we will still introduce several available French and EU crop databases, which could easily be treated in the same way as the grain maize application.

Chapter 2

Analysis of available datasets

Chapter 2 Analysis of available datasets

The second database is the weather and climate data. For the current or historical record, we will describe the reanalysis data and introduce some climate variables as well as their pre-processing. Then, the climate simulations, which provide the future climate data, are described. Finally, other available climate datasets (i.e. observations) are discussed.

Agricultural data

Coffee database

Coffee is one of the most traded agricultural commodities in the world. It is a major cash crop and an important foreign exchange earner for many coffeeproducing countries. The two most wildly cultivated coffee species are Arabica (Coffea arabica) and Robusta (Coffea canephora Pierre ex A. Froehner).

The global coffee cultivation

Coffee is cultivated in over 70 countries ranging in latitude from 30 • N to 30 • S. Arabica producing areas require a climate with average annual temperatures between 18 and 21 • C, a rainfall pattern between 1100 and 2000 mm [START_REF] Damatta | Impacts of drought and temperature stress on coffee[END_REF], over elevations typically from 950 to 1950 m [START_REF] Davis | An annotated taxonomic conspectus of the genus Coffea (Rubiaceae)[END_REF]. Compared to Arabica, Robusta is successfully grown at the elevation between 50 and 1500 m [START_REF] Davis | An annotated taxonomic conspectus of the genus Coffea (Rubiaceae)[END_REF], with annual temperatures from 22 • C to 26 • C [ [START_REF] Damatta | Impacts of drought and temperature stress on coffee[END_REF], and with an annual rainfall pattern between 1200 and 2500 mm.

Figure 2.1 shows the global map of major coffee-growing regions. The top 11 coffee-producing countries are also given; they account for 85 % of the global coffee production in 2019 [START_REF] Usda | Coffee: World markets and trade[END_REF]FAO, 2019]. Brazil is the leading coffee producer, with over one-third of the world total in 2019 (about 43 % of Arabica and 20 % of Robusta coffee) [START_REF] Usda | Coffee: World markets and trade[END_REF]FAO, 2019;of Agriculture Foreign Agricultural Service, 2018]. It is followed by Vietnam with about 20 % of total production and more than 40 % of Robusta coffee [START_REF] Usda | Coffee: World markets and trade[END_REF]FAO, 2019]. Since most of the coffee-producing countries are located around the Equator, they have pretty similar crop years (i.e. the 12-month period from one harvest to the next): from April to March (e.g. Brazil, Indonesia, and Peru), from October to September (e.g. Vietnam, Ethiopia, and India) [START_REF] Marsh | Diversification by smallholder farmers: Viet Nam Robusta Coffee[END_REF]. However, many countries have more than one crop year (several coffee harvest seasons). For instance, in Colombia, besides the main harvest (October -December), there is a second harvest (April -July), especially in the southeast regions of Colombia (see more at https:// federaciondecafeteros.org/wp/cosecha-cafetera/). This makes the modelling of coffee yield even more difficult.

No Country

ID Main crop year [START_REF] Marsh | Diversification by smallholder farmers: Viet Nam Robusta Coffee[END_REF] Coffee type [START_REF] Marsh | Diversification by smallholder farmers: Viet Nam Robusta Coffee[END_REF] Coffee Production (Metric Tons) Chapter 2 Analysis of available datasets

Administrative division

Let us first introduce the administrative division of some major coffee-producing countries. The administrative division [www.cia.gov., 2021] or administrative unit [FAO, 2015] are generic names for geographical areas into which a country is divided, and they can be named differently. The common names for the principal (or largest) administrative divisions (or sub-national level) are states, provinces, departments, regions, cantons, etc. These principal units are often subdivided into smaller ones known by names such as circuits, counties, municipalities, districts, etc. And these division units are named quite differently among coffee-producing countries. For instance, Brazil is divided into 26 states (plus 1 federal district) and about 5568 municipalities (plus 1 state district), while Vietnam consists of 58 provinces (plus 5 municipalities) and 585 districts (plus several provincial cities and islands). To simplify these units, we will use the terms: the first administrative level L1 for the largest unit and the second administrative level L2 for the next subdivision. In the following, these two division levels will be used for all the coffee-related applications. Smaller divisions (e.g. third or fourth levels) will not be considered in this thesis for practical reasons (such as the resolution of climate models).

Available coffee databases worldwide

Available datasets at the global scale only provided information on coffee production at the country level, such as in the Food and Agriculture Organisation (FAO), the US Department of Agriculture (USDA), or the International Coffee Organisation (ICO). To our knowledge, there is no existing global dataset of coffee-producing regions at the sub-national levels (e.g. L1 and L2) for coffee yield (or production, area). Thus, we gathered as much coffee production information as possible, from various sources, mainly from the Ministry of Agriculture and the National Statistics Office of the major coffee-producing countries. Gathering this database is challenging for many reasons: the coffee data are of a commercial nature, the (freely available) data are in various formats and languages, and we need to contact several agencies. The database we have collected is a very nice resource now that could be shared if needed for collaborative projects. Table 2.1 summarises the collected data on production, harvested and/or cultivated area, yield, and price. The timeline of available yield datasets is also presented (see Tab. 2.1). We gathered data for several major coffee-growing countries such as Brazil, Vietnam, Colombia, Indonesia, and India. The coffee data of the two leading coffee producers (Brazil and Vietnam) are available over longer records (still not long enough) than in other countries at different administrative levels. We also obtained long temporal records (25 years of yield data) for coffee in Peru, but these data are only available at the first administrative level. Other datasets with less than ten years (e.g. the case of Indonesia, India, Guatemala, and Tanzania) seem insufficient to study the link with weather directly. total coffee production, followed by Espirito Santo, Sao Paulo, Bahia, and Rondonia, with about 20.05 %, 10.50 %, 8.03 %, and 4 % of the total production of Brazil (IBGE). These five states represent more than 93 % of Brazil's production. Although these states cultivate both Robusta and Arabica, Arabica is predominant. For example, only Arabica coffee is cultivated in Minas Gerais, contributing to about half of the country's Arabica production. In Brazil, most coffee plantations are harvested in the dry seasons of June through September [START_REF] Souza | Plant-Parasitic Nematodes of Coffee[END_REF].

For Vietnam, the leading Robusta coffee producer, most of its coffee is grown in the Central Highlands, which includes five provinces: Dak Lak, Lam Dong, Dak Nong, Gia Lai, and Kon Tum (shown in Fig. 2.3b). From 2000 to 2018, the Central Highlands produced 93.01 % of Vietnam's total coffee production: on average, 35.36 % for Dak Lak, 27.33 % for Lam Dong, 14.82 % for Dak Nong, 13.27 % for Gia Lai, and 2.23 % for Kon Tum. Central Highlands lies on a series of contiguous plateaus with an elevation ranging from 500 to 1500 m. The coffee harvest season occurs from October to January. The coffee yield time series of Brazil and Vietnam, for both L1 and L2 levels, are shown in Fig. 2.4. In detail, Fig. 2.4a1 presents the yield over five major coffee-producing states in Brazil. Generally, the yields show an increasing trend from the 1990s. They suddenly drop in 2002 and then start to increase again. This drop could be due to the drop in Brazilian Arabica price at the end of the twentieth-century [ICO, 2019]. Due to low prices, the farmers did not apply or reduce management practices (e.g. no pruning, less fertiliser), resulting in low coffee yields in the following years. Figure 2.4a2 describes the yield data of 100 municipalities in Minas Gerais (Brazil). These datasets (Brazilian coffee data at the second administrative level in general) are very discrete and not well organised: yield data in many municipalities are only available for some short periods and not continuous. In addition, Fig. 2.4b1 presents the coffee yield over five major coffee-producing provinces in Vietnam, and Fig. 2.4b2 for several districts in Dak Lak. Overall, the coffee yield data at both L1 and L2 levels steadily increased from 2000 to 2018.

Other coffee-related data

Together with the coffee information and climate data (which will be introduced in Sect. 2.2), we also need other datasets such as elevation and soil type. These two factors have been proved to impact the coffee quantity and quality [START_REF] Tsegay | Effect of altitude of coffee plants on the composition of fatty acids of green coffee beans[END_REF][START_REF] Getachew | The relationship between elevation, soil temperatures, soil chemical characteristics, and green coffee bean quality and biochemistry in southwest ethiopia[END_REF].

Elevation

Digital elevation data, SRTM 1 Arc-Second Global [USGS EROS Center, 2018] data, are downloaded from https://www.usgs.gov/centers/eros/science/ usgs-eros-archive-digital-elevation-shuttle-radar-topographymission-srtm-1. This dataset offers worldwide coverage of void filled data at a resolution of 1 arc-second (30 m) and provides open distribution of this highresolution global data set. Figure 2.5 presents the elevation data over Brazil and Vietnam.

Soil type

The harmonised world soil database [START_REF] Fischer | Global agro-ecological zones assessment for agriculture[END_REF] is the result of a collaboration between the FAO with IIASA, ISRIC-World Soil Information, Institute of Soil Science, Chinese Academy of Sciences (ISSCAS), and the Joint Research Centre of the European Commission (JRC). It provides a 30 arc-second raster database with over 15 000 different soil mapping units. The data is available for download at https://www.fao.org/soils-portal/soil-survey/soilmaps-and-databases/harmonized-world-soil-database-v12/ru/. Examples of soil data over Brazil and Vietnam regions are shown in Fig. 2.6.

European Union (French) crops databases

The NUTS classification

In addition to the classical administrative division (as presented in Sect. 2.1.1.2), the nomenclature of territorial units for statistics, abbreviated as NUTS (from the French "nomenclature des unités territoriales statistiques"), is also commonly In the following of this thesis, we will use this NUTS classification for the study related to the EU area. For instance, the crop database in Sect. 2.1.2.4 follows the NUTS version 2016, which are detailed in https://eur-lex.europa.eu/ legal-content/EN/TXT/HTML/?uri=CELEX:32016R2066&from=FR.

General information of the French agriculture

Agriculture in France 1 benefits from a large agricultural area (occupying about 53.2 % of the land surface) and a favourable geographical and climatic situation. The main products are cereals (wheat, maize), sugar, wine, milk, dairy products, fruit and vegetables, livestock, and meat. According to the FAO (data in 2018), France is among the leading countries in many agricultural products, for instance, the 5 th world producer of wheat, the 2 nd world producer of mixed cereals, the 6 th world producer of whole cereals, the 8 th world producer of fresh corn, the 2 nd world producer of barley, the 5 th world producer of grapes, etc.

The administrative divisions of France are organized into 13 regions, including Corsica [Ministère de la Cohésion des territoires et des Relations avec les collectivités territoriales, 2021] in 2021. The regions are subdivided into 96 departments (in French: département), which corresponds to the third level of NUTS (i.e. NUTS 3). The borders of French departments have changed over time [START_REF] Schauberger | Yield trends, variability and stagnation analysis of major crops in France over more than a century[END_REF][START_REF] Ceglar | Time-varying impact of climate on maize and wheat yields in France since[END_REF]. In this thesis, data are gathered for 96 French mainland departments in their current form (shown in Fig. 2.7), and the historical yield values were subsumed to modern departments.

Available French crops databases

The French crops data (area, production, and yield) on the regional level (department) were collected from the Agreste website (https://agreste.agriculture.gouv.fr, last access: 15 December 2020) compiled by the Service de la Statistique et de la Prospective (SSP), Secrétariat Général du Ministère de l'Agriculture, de l'Agroalimentaire et de la Forêt (MAAF), France). The available data are presented in Tab. 2.2. For example, Fig. 2.8 shows yield time series of 96 French departments from 1989 to 2010 for different crops (winter wheat, winter durum wheat, winter oat, and spring oat).

In addition, we also found long time series data for several crops in France in the literature (see Tab. 2.3). These data were gathered by [START_REF] Schauberger | Yield trends, variability and stagnation analysis of major crops in France over more than a century[END_REF] from several sources, for instance, books of national agricultural statistics ("Statistique agricole annuelle") compiled by the French Ministry of Agriculture, digital statistics from the Agreste website.

Available EU crops databases

Like coffee, several services provide global data at the country level for many crops, such as FAO (https://www.fao.org/faostat/en/#data/QCL), USDA (https://www.usda.gov/topics/data). For the EU countries, there is the European statistical office or Eurostat (https://ec.europa.eu/eurostat/ data/database), which offers a range of important and high-quality data. However, it often provides data only at the country level. At higher levels or subnational levels, the data availability depends on the types of crops or countries.

Among many types of crops produced, wheat is one of the most important crops in EU countries. The EU produced about 255 million metric tons of wheat • total wheat (Eurostat code C1100),

• soft wheat (C1110),

• durum wheat (C1120).

Figure 2.9 shows the total wheat yield data (in t•ha -1 ): (a) yield in 2015 over 27 EU countries, and (b) yield time series of 16 German regions (i.e. federated states or Bundesländer in German) at NUTS 1 level.

For a long time, countries have tried to collect and organize data to have national agricultural statistics, and institutions have attempted to obtain data on a global scale. However, data collection is always a challenging task: there is no data in some regions and/or some periods, the data are not uniform, or data were not stored in the past, etc. We need to accept all the limitations and thus take them into account in how we will use the available and accessible database.

Weather and climate data

As our objective is to study weather and climate impacts on agriculture, current (or historical) and future climate data will be required.

Historical record

What is climate reanalysis?

Weather records for the present are quite complete, with many different modern observation methods. These observation data, however, are not evenly distributed around the globe. Data availability for the past is much more scarce. The reanalysis, a scientific method for developing a comprehensive climate record, aims at solving this problem. Climate reanalysis combines past observations with climate models to deliver a complete and consistent picture of the past climate as close to reality as possible, extending back several decades or more. It contains estimates of atmospheric parameters (e.g. air temperature, pressure, and wind at different altitudes) and surface parameters (e.g. rainfall, soil moisture content, ocean-wave height, and sea-surface temperature). The estimates are produced globally, and they span an extended period.

Various climate reanalysis products are freely available, for instance, the NCEP/N-CAR reanalysis [START_REF] Kalnay | The ncep/ncar 40-year reanalysis project[END_REF] (a joint product from the National Centers for Environmental Prediction and the National Center for Atmospheric Research), the MERRA reanalysis [START_REF] Gelaro | The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2)[END_REF] (Modern-Era Retrospective analysis for Research and Applications), and the ECMWF reanalysis [START_REF] Hersbach | Operational global reanalysis: progress, future directions and synergies with nwp[END_REF] (European Center for Medium-Range Weather Forecasts).

Reanalysis products are extensively used in climate research and applications. Among many available products, many studies have shown that the ECMWF reanalysis is suitable for crop yield simulations [START_REF] Challinor | Simulation of crop yields using ERA-40: Limits to skill and nonstationarity in weather-yield relationships[END_REF][START_REF] Mathieu | Modèles d'impact statistiques en agriculture : de la prévision saisonnière à la prévision à long terme, en passant par les estimations annuelles[END_REF]. This thesis uses ERA5-Land covering the land areas of the ERA5 reanalysis [START_REF] Hersbach | Operational global reanalysis: progress, future directions and synergies with nwp[END_REF], which is the ECMWF latest2 climate reanalysis from 1950 onward. The native spatial resolution of this dataset is 9 km on a reduced Gaussian grid (TCo1279). Then the data are regridded to a regular lat-lon grid of 0.1 • × 0.1 • (about 11 km × 11 km at the Equator). These 0.1 • × 0.1 • gridded data are used in this thesis and can be obtained from https: //cds.climate.copernicus.eu/cdsapp#!/search?type=dataset.

Climate variables and their pre-processing

The monthly means of different weather variables were collected from the ERA5-Land, e.g. total precipitation, 2 m mean temperature, solar radiation, and evaporation. Figure 2.10 presents the annual total precipitation (in 10 3 mm) and the average 2 m temperature (in • C) at a global scale, in 2000.

In this thesis, some applications can use the gridded climate data directly. However, some applications may require matching the climate resolution to the one of the crop yield data, i.e. in the yield estimation tasks (see in chapter 4). In these cases, the gridded climate data are projected from its original 0.1 • × 0.1 • grid into the crop administrative levels (for example, the provincial and/or district levels for Robusta coffee in Vietnam (e.g. Fig 2.11), the department level for crops in France, and different NUTS levels for crops over Europe). This has been done by aggregating the gridded data over the administrative shapes as follows:

(1) if the shape is smaller than the cell, the gridded value will be representative of the region; (2) if the shape includes several cells, the weather data will be averaged based on the area of cells inside the shape. In addition, climate data are used to derive several agro-climatic variables, which will be described later in the corresponding chapters (e.g. in chapter 7).

Climate simulations

Climate models are numerical representations of the climate system based on the fundamental laws of physics, fluid dynamics, and thermodynamics. These models characterise how energy and matter interact in different parts of the earth (e.g. atmosphere, ocean, land). Climate simulations are obtained by running these numerical models, which may cover either the entire globe (i.e. global climate models (GCMs) or also known as general circulation models) or a specific region (i.e. regional climate models (RCMs)). Climate projections are the simulations in future decades (typically until 2100) based on assumed emissions scenarios, i.e. scenarios for the concentrations of greenhouse gases, aerosols, and other atmospheric constituents that affect the planet's radiative balance. These simulations are among our primary means to understand how the climate has changed in the past and may change in the future. In the following, we will first describe several future scenarios that drive the GCMs and RCMs. Then, these two types of models, as well as their limitations, will be introduced.

Future scenarios

Over the past few years, a new set of scenarios (i.e. pathways) has been developed by a joint community effort (including climate scientists, economists, and energy system modellers). These scenarios are alternative possibilities that could happen in the future, showing the relationship between human activities, emissions, concentrations, and climate change. The First Intergovernmental Panel on Climate Change Assessment Report (IPCC FAR) in 1990 discussed four emissions scenarios. Today, these scenarios have expanded to encompass various timedependent or transient scenarios and adapt to some significant changes (e.g. populations, economy) in the past 20 years. The most recent set of time-dependent scenarios is the "Representative Concentration Pathways" (RCPs) [START_REF] Moss | Towards New Scenarios for Analysis of Emissions, Climate Change, Impacts, and Response Strategies[END_REF][START_REF] Moss | The next generation of scenarios for climate change research and assessment[END_REF], describing different levels of greenhouse gases and other radiative forcings that might occur in the future. RCPs include four pathways which are selected and defined by their total radiative forcing in 2100 (2.6, 4.5, 6.0, and 8.5 in W•m -2 ) [START_REF] Van Vuuren | The representative concentration pathways: an overview[END_REF]. Figure 2.12, adopted from the study of [START_REF] Meinshausen | The rcp greenhouse gas concentrations and their extensions from 1765 to 2300[END_REF], presents harmonised emissions under the four RCP scenarios.

• The RCP2.6 scenario or a stringent mitigation scenario is a so-called "peak" scenario, which implies that the radiative forcing level reaches 3.1 W•m -2 by mid-century but returns to 2.6 W•m -2 at the end of the century.

• The RCP4.5 and RCP6.0 scenarios are intermediate scenarios. Emissions in RCP4.5 peak around 2040 and in RCP6.0 peak around 2080, then decline (Fig. 2.12).

• The RCP8.5 is a high greenhouse gas concentration scenario in which radiative forcing reaches greater than 8.5 W•m -2 by 2100 and continues to rise for some amount of time. It is generally considered as the basis for worst-case climate change scenarios.

(

( ) ( ) ( ) ( ) )
FIGURE 2.12: Harmonised emissions under the four RCP scenarios. The nonharmonised scenarios (dashed lines) are in most cases marginally different from the harmonised emissions (solid lines). For illustrative purposes, emissions are weighted with IPCC SAR global warming potentials of a 100-year time horizon (IPCC 1996). The figure is adopted from Fig. 2 in [START_REF] Meinshausen | The rcp greenhouse gas concentrations and their extensions from 1765 to 2300[END_REF].

Although these four RCPs are consistent with certain socioeconomic assumptions, they did not include any socioeconomic data. A new set of scenarios was thus developed, the "Shared Socioeconomic Pathways" (SSPs), including information on population, economic growth, education, urbanisation, and the rate of technological development [START_REF] O'neill | The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6[END_REF][START_REF] Riahi | The shared socioeconomic pathways and their energy, land use, and greenhouse gas emissions implications: An overview[END_REF]. These SSPs are based on five narratives describing alternative socioeconomic developments [START_REF] O'neill | A new scenario framework for climate change research: the concept of shared socioeconomic pathways[END_REF][START_REF] O'neill | The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6[END_REF][START_REF] Riahi | The shared socioeconomic pathways and their energy, land use, and greenhouse gas emissions implications: An overview[END_REF], which are:

• SSP1: sustainability -taking the green road (low challenges to mitigation and adaptation),

• SSP2: middle of the road (medium challenges to mitigation and adaptation),

• SSP3: regional rivalry -a rocky road (high challenges to mitigation and adaptation),

• SSP4: inequality -a road divided (low challenges to mitigation, high challenges to adaptation),

• SSP5: fossil-fuelled development -taking the highway (high challenges to mitigation, low challenges to adaptation).

Global climate models (GCMs)

The Coupled Model Intercomparison Project (CMIP) is a project of the World Climate Research Programme (WCRP)'s Working Group of Coupled Modelling (WGCM). This project started 20 years ago as a comparison of a handful of early global coupled climate models involving multiple international modelling teams worldwide [START_REF] Meehl | Intercomparison makes for a better climate model[END_REF][START_REF] Eyring | Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization[END_REF]. The objective of CMIP is to better understand past, present, and future climate change and variability in a multi-model framework. CMIP has developed in phases, with the simulations of the fifth phase (CMIP5) completed and the sixth phase (CMIP6) about to accomplish3 . The main difference between these two phases is that CMIP5 projections are available based on RCP scenarios, while CMIP6 uses socioeconomic pathways (SSPs). Also, CMIP6 is expected to produce considerably more data than CMIP5 due to several factors: the increases in model resolution, an increase in the number of participating modelling centres, and an increased number of experiments from 97 in CMIP5 to 321 experiments in CMIP6 [START_REF] Balaji | Requirements for a global data infrastructure in support of cmip6[END_REF][START_REF] Petrie | Coordinating an operational data distribution network for cmip6 data[END_REF].

In the following, we will use CMIP6 data as the inputs to review different bias correction methods (i.e. calibration methods) of climate models (chapter 6). In detail, the monthly mean of temperature and precipitation data from CNRM-CM6-1-HR [START_REF] Voldoire | Cnrm-cerfacs cnrm-cm6-1-hr model output prepared for cmip6 highresmip[END_REF][START_REF] Voldoire | Evaluation of cmip6 deck experiments with cnrm-cm6-1[END_REF] are considered. CNRM-CM6-1-HR is the high resolution of the fully coupled atmosphere-ocean general circulation model of sixth-generation jointly developed by CNRM (Centre National de Recherches Météorologiques) and Cerfacs (Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique). These data belong to the highresolution model intercomparison project of the CMIP6 [START_REF] Eyring | Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization[END_REF] and can be downloaded from https://esgf-node.ipsl.upmc.fr/projects/ cmip6-ipsl/. They are available from 1850 to 2100 at a resolution of 1 • × 1 • . Due to time limitations, we consider here only one model. Other climate data (e.g. from different institutes such as IPSL, GERICS, or MPI) or even an ensemble of models can also be considered, for example, in future works.

Regional climate models (RCMs)

In practice, a climate impact model (e.g. on agriculture) often requires more detailed climate data than the ones offered by GCMs. Moreover, GCMs could introduce large uncertainties, especially over some specific regions with vastly differing landscapes, e.g. Southeast Asia -one of the world's most vulnerable regions to climate change. Thus, regional climate models (RCMs), applied over a limited area and driven by GCMs, could actually be considered as a better alternative source for climate change impact studies [START_REF] Nguyen | More intense daily precipitation in cordex-sea regional climate models than their forcing global climate models over southeast asia[END_REF].

In chapter 7, we will investigate the climate change impacts on coffee production using the RCMs from the WCRP's Coordinated Regional Climate Downscaling Experiment (CORDEX). CORDEX initiative delivered downscaled simulations for various GCMs from CMIP5 [START_REF] Meehl | Decadal prediction: Can it be skillful[END_REF][START_REF] Taylor | An overview of cmip5 and the experiment design[END_REF] to higher resolution regional models for 14 regions worldwide. Here, we will focus on two regions covering two leading coffee-producing countries: South America (SAM) and Southeast Asia (SEA). Figure 2.13 presents detailed descriptions of these two regions. (4) number of point in the North-South direction: Ny.

The regional climate data are obtained from several model simulations of CORDEX, with 0.22 • × 0.22 • resolution, based on RCP scenarios. The data are available at https://esg-dn1.nsc.liu.se/search/cordex/. Details of models used will be introduced later in chapter 7.

Climate simulations are uncertain

Climate simulations provide information on future climate at regional and global scales. They are our unique source of knowledge about climate change and thus are widely used in climate impact models. However, they often suffer from substantial biases and errors compared to the real-world climate described by observations [START_REF] Flato | Evaluation of climate models[END_REF][START_REF] Kotlarski | Regional climate modeling on european scales: a joint standard evaluation of the eurocordex rcm ensemble[END_REF]. This imperfection often comes from two main types of uncertainty: structural and parametric. First, structural uncertainty occurs because models do not represent all essential physical processes correctly. Second, errors in parameterization cause parametric uncertainty. This is a method that replaces processes that occur at too small-scale (spatial or temporal) to be physically represented in the model by a simplified process. In the climate model, these parameterizations are roughly accurate over the scale of a grid box (e.g. 100 × 100 km). In addition, different models are based on different patterns of natural variability, and this variability is chaotic. Last but not least, future scenarios are driven by human activities which are impossible to predict. The feedbacks (e.g. cloud feedback) add uncertainty to climate models.

It is generally difficult to identify one type of uncertainty that matters most. It rather depends on the study period. For instance, for short-term projections (e.g. the next 10-20 years), the natural internal variability is the most important source of uncertainty [Hawkins andSutton, 2009, 2011]. Thus, it is better to consider as many climate model simulations as possible, to reduce/limit the total uncertainty. On the other hand, for beyond 20-30 years, human choices (i.e. different model responses to prescribed emissions) are most important for temperaturerelated impacts and for drought and extreme precipitation over larger spatial scales [Hawkins andSutton, 2009, 2011]. Therefore, it is better to use a range of future scenarios that encompass possible outcomes in this case.

Other datasets

In addition to the reanalysis data, in situ data or station data exists over our study area, for instance, in the coffee-growing area in Vietnam. The impact models will probably be more precise with better weather data, i.e. the in situ observations. Nevertheless, this thesis will consider the reanalysis data (ERA5-Land) as historical records for several reasons.

Firstly, although station data exist, they are often not open access. For example, there are at least 13 meteorological stations (shown in Fig. 2.14) that are installed and in operational mode across the Central Highlands of Vietnam [START_REF] Ngo-Thanh | A distinction between summer rainy season and summer monsoon season over the central highlands of vietnam[END_REF][START_REF] Pham-Thanh | Predictability of the rainy season onset date in central highlands of vietnam[END_REF]. However, among these 13 stations, we can (freely) access the data of only two stations via the WMO GTS (World Meteorological Organization's Global Telecommunication System): Kon Tum and Pleiku.

The availability of in situ measurements can be improved in other regions. However, the most important reason is that we aim to apply the impact models for not only current climate but also climate change experiments. "Essentially, all models are wrong, but some are useful." [START_REF] Box | Robustness in the strategy of scientific model building[END_REF] As introduced in chapter 1, this thesis aims at measuring the crop's sensitivity to weather and/or climate conditions by developing statistical crop yield models based on weather information. In order to accomplish this objective, we need to identify the appropriate predictors, the adequate statistical models, and estimate their true generalization quality. In this chapter, we will introduce basic concepts and various types of commonly used statistical and machine learning models. Then, the model complexity and overfitting problem are discussed. Several metrics used to evaluate the model performance are also presented. After that, we will present several techniques to avoid overfitting, e.g. regularization techniques and dimensionality reduction. Finally, the cross-validation techniques, including a newly proposed "leave-two-out" nested cross-validation method, will be discussed, especially in the context of having limited data. This proposed technique is the subject of a publication in Geoscientific Model Development journal entitled "Nested leave-two-out cross-validation for the optimal crop yield model selection" [START_REF] Dinh | Climate change impacts on robusta coffee production in vietnam[END_REF].

Although it is a methodology chapter, we rely here on the example of yield modelling or coffee suitability classification to facilitate understanding and begin to introduce the problems that will be studied in the following chapters.

Statistical/machine learning models

Statistical or machine learning models can broadly be categorized as unsupervised or supervised.

Unsupervised learnings refer to the use of machine learning algorithms to identify patterns in datasets containing data points that are neither classified nor labelled. Two main methods used in unsupervised learning include clustering and dimensionality reduction.

In contrast, supervised learnings involve learning a function that maps an input to an output based on example input-output pairs [START_REF] Russell | Artificial intelligence : a modern approach[END_REF]. Regression and classification are the two main supervised learning algorithms. Both of these algorithms are used for prediction and work with the labelled datasets. The main difference between these two approaches is how they are used for different applications. Regression algorithms are used to predict continuous values (e.g. crop yield). Meanwhile, classification algorithms are used to predict or classify discrete values (e.g. suitable (class 1) versus unsuitable (class 0)). In practice, many types of algorithms or models can be used to accomplish these two tasks, for instance:

• regression: linear regression, decision tree, random forest, neural networks, mixed-effects model, etc.;

• classification: logistic regression, decision tree, random forest, support vector machines, Naïve Bayes, neural networks, etc.

This thesis focuses on supervised learning. Supervised learning will be used in the following chapters for several applications, for example, predicting the crop yield or classifying the suitable crop areas. To do these tasks, the main idea is to use an algorithm to learn the mapping function from the input (i.e. one or more explanatory variables or predictors X) to the output (i.e. a dependent variable or a quantitative response Y ), which can be noted as:

Y = f (X) + (3.1)
where f is some fixed but unknown function of X and is a random error term. The function f can be based on several models depending on the complexity of the application. In the following, we will present several common algorithms such as linear regression, random forest, neural networks, and mixedeffects models. Some of them will be applied to real-world applications presented in chapters 4, 5, and 7.

Linear regression

The model The linear regression (LIN) model is known as the most simple model and it is probably the most used model among supervised learning models. It is based on a linear relationship between a scalar response (dependent variable) and one or more explanatory variables (or also called predictors in following chapters). In our study, we will consider several explanatory variables, the model is then known as multiple linear regression. The relationship between the crop yield data1 y t ( t = 1, 2, • • • , n samp corresponds to the year) and the meteorological variables x tj (j = 1, 2, • • • , n input where n input is the number of input variables2 ) is formulated as:

y t = β 0 + β 1 • x t1 + • • • + β p • x tn input + t , (3.2)
where β i is the regression coefficients and is a disturbance term or error variable. This represents the errors of the model, typically by using a random variable. These n samp equations, corresponding to the n samp statistical samples3 , can be written in matrix notation as:

Y = Xβ + (3.3) where Y =       y 1 y 2 . . . y nsamp       ∈ R nsamp , X =       1 x 1,1 . . . x 1,n input 1 x 2,1 . . . x 2,n input . . . . . . . . . . . . 1 x n,1 . . . x nsamp,n input       ∈ R nsamp×(n input +1) , β =         β 0 β 1 β 2 . . . β n input         ∈ R n input +1 , =       1 2 . . . nsamp       ∈ R nsamp .
is a random vector that follows a multivariate normal distribution, ∼ N (0, σ 2 0 I n ).

Estimation method

The most common procedure for parameter estimation in linear regression is the least-squares estimation. In this least-squares setting, the optimum parameters are defined as such that they minimise the sum of mean squared loss or loss function or the residual sum of squares (RSS):

β = arg min β RSS = arg min β ||Y - Xβ|| 2
The loss function can be rewritten as:

RSS = ||Y -Xβ|| 2 = (Y -Xβ) T (Y -Xβ) = Y T Y -Y T Xβ -β T X T Y + β T X T Xβ (3.4)
The gradient of the loss function with respect to the parameters β is:

∂RSS ∂β = ∂(Y T Y -Y T Xβ -β T X T Y + β T X T Xβ) ∂β = -2X T Y + 2X T Xβ (3.5)
The optimum parameter is obtained when the gradient of the loss function is set to zero:

-2X T Y + 2X T Xβ = 0 =⇒ X T Y = X T Xβ =⇒ β = (X T X) -1 X T Y (3.6)
Note that ((X T X) -1 X T ) is the pseudo-inverse of the rectangular matrix X.

Random forest

Random forest (RF) [START_REF] Breiman | Random forests[END_REF] is a non-parametric machine learning algorithm which can be applied for regression (e.g. predict crop yield in [START_REF] Jeong | Random forests for global and regional crop yield predictions[END_REF]; [START_REF] Beillouin | Impact of extreme weather conditions on european crop production in 2018. Philosophical transactions of the Royal Society of London[END_REF]) and classification (e.g. crop classification in [START_REF] Tatsumi | Crop classification of upland fields using random forest of time-series landsat 7 etm+ data[END_REF]). The "forest" it builds is an ensemble of decision trees, usually trained with the "bagging" method [START_REF] Breiman | Bagging predictors[END_REF]. The general idea of the bagging method is that a combination of learning models increases the overall result. A RF algorithm has several tuning parameters that need to be set before training, in particular the node size, the number of trees, and the number of features.

In general, when using RF for regression, the forest predicts by taking the average or mean of the outputs from all trees. On the other hand, when using RF for classification, each tree gives a "vote". Then, the forest chooses the classification with the majority of the "votes".

Neural networks

Neural networks (NN)4 models are non-linear statistical models [START_REF] Haykin | Neural Networks: A Comprehensive Foundation[END_REF][START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF]. Its name is inspired by the biological neural networks that structure animal brains. In practice, these models are extremely efficient but they often require a sufficient amount of data.

Artificial neurons

Like a biological neuron (Fig. 3.1a) which has synapses (connection with other neurons), dendrites ("neuron input"), an axon ("neuron output"), and a nucleus (activating the output according to input stimulation), an artificial neuron (Fig. 3.1b) is a model characterised by: (1) the inputs x 1 , • • • , x n , (2) a vector of weighs W , and (3) an activation function f : R → R, and an output y.

The activation function f applies to the weighted sum of all the inputs, whose coefficients are a weight vector W = [w 0 , • • • w n ] ∈ R n+1 associated with each input. We call b the bias (or "offset", or "threshold") of the neuron. The values of the weights are estimated during the training. The output y can be finally presented as:

y = f n i=0 w i x i + b .
(3.7)

There are numerous possible activation functions:

• the linear function: f is identity;

• the simple logistic/sigmoid function: • the rectified linear unit (or ReLU) function: f (x) = max(0, x);

f (x) = 1 1+e -x ;
• the tanh function (hyperbolic tangent):

f (x) = (e x -e -x ) (e x +e -x )
In practice, the used activation function are often continuous and differentiable with respect to the weights W for optimisation purposes. The most commonly used activation function in a NN is the sigmoid function.

Neural networks architecture

The neurons are aggregated into layers. Neurons of one layer connect only to neurons of the immediately preceding and immediately following layers. These layers form different parts of a NN, typically three parts:

• an input layer: receives external data (i.e. inputs);

• one or more hidden layers: is the intermediate layers where neurons take in a set of weighted inputs and produce an output though an activation function;

• an output layer: produces the ultimate results performed by the network to the outer world. The simplest networks contain no hidden layers and are equivalent to linear regression.

Types of NN and feedforward neural networks

Neural networks can be categorized into different types based on their structure, data flow, neurons used and their density, layers and their depth activation filters, etc. Some common types of NN include:

• feedforward neural networks;

• perceptron and multilayer perceptron neural networks;

• radial basis function neural networks;

• recurrent neural networks.

In this thesis, we will study one of the simplest types of NN which is feedforward neural networks [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF][START_REF] Schmidhuber | Deep learning in neural networks: An overview[END_REF]. In this network, there is only one direction forward from the input nodes: through the hidden nodes and to the output nodes. Only one hidden layer with n neuron neurons is considered in the architecture here (Fig. 3.2). The outputs of the nodes in one layer are inputs to the next layer. The inputs to each node are combined using a weighted linear combination. It is then transferred to an activation function of the nonlinear sigmoid function σ. The crop yield data y is modeled in this case by the following equation:

y = nneuron j=1 w j × σ n input p=1 w jp x p + b hidden + b output (3.8)
where w represents the weights associated to each connection in the NN architecture, the x i 's are the meteorological variables, and b are the NN biases.

Linear mixed-effects

While classic linear models handle between-subjects (or between-groups) study designs, mixed-effects (ME) models5 can handle both between-subjects and withinsubjects study designs. In other words, ME models are often used in settings where repeated measurements are made of the same statistical units over time (i.e. longitudinal data) or where measurements are made on clusters of related statistical units (i.e. grouped data).

A ME model contains both "fixed" and "random" effects. The term "effect" is related to the model's response to an input perturbation. Fixed effects are variables that are constant across individuals; these variables do not change or change at a constant rate over time. Estimating these effects is the traditional domain of regression modelling. Random effects, on the other hand, are sample-dependent random variables; these variables are random and unpredictable. These random terms are useful when data fall into natural groups, e.g. soil types, geographic categories, or land-use groups [START_REF] Mathieu | Statistical weather-impact models: An application of neural networks and mixed effects for corn production over the United States[END_REF]. In the case of studying agricultural yield determinants, both fixed and random effects are typically present. For instance, the practices (types of fertilisers or herbicides) are fixed, but geographic locality or years may be random.

Adding random effects to a model extends the reliability of the inferences by considering the variability between groups. Taking into account the group information, a ME model thus can make a compromise between the complete pooling and no pooling methods. The complete pooling method means that the model ignores variations between groups and treats all observations as part of the same group or pool, resulting in a general model. While no pooling method assumes groups share no similarities and fits each group with a separate model.

Mixed-effects models are applied in many disciplines where multiple correlated measurements are made on each unit of interest. We first start with some simple linear mixed-effects models [START_REF] Lindstrom | Newton-raphson and em algorithms for linear mixed-effects models for repeated-measures data[END_REF]]. Let's assume that the (yield) data, used to calibrate an impact model, are grouped into m distinct groups, e.g. m agro-ecological zones. We note here y it the t th observation of group

i (t = 1, 2, • • • , n i ; i = 1, 2, • • • , m
) that we are trying to predict, and x it is the t th value in group i of a predictor x (only one predictor is considered in this first simple approach). A linear mixed-effects model can be written as:

y it = β 0i + β 1 x it + it , (3.9)
where β 0i presents an intercept for each group, and it can be rewritten as:

β 0i = β 00 + b 0i
, where b 0i is assumed to be normally distributed b 0i ∼ N (0, σ 2 0 ). The term it is the measurement or process error: it ∼ N (0, σ 2 )). b 0i and it are independent to each other. In a real application, y it can be, for instance, the coffee yield anomaly, t is the year from 2000 to 2019 (t = 2000, 2001, • • • , 2019), and i = 1, 2, • • • , 10) represent 10 agro-ecological zones/groups.

By reordering the parameters depending on the group, we obtain:

y it = β 00 + β 1 x it fixed effects + b 0i random effects + it . (3.10)
The parameter related to the explanatory variable x (i.e. slope) can also depend on the group and thus have a random effect:

y it = β 0i + β 1i x it + it .
(3.11)

Again, we can rewrite β 0i and β 1i as:

β 0i = β 00 + b 0i , b 0i ∼ N (0, σ 2 0 ), β 1i = β 10 + b 1i , b 1i ∼ N (0, σ 2 1 ), or generally: b i = b 0i b 1i ∼ N 0, σ 2 0 0 0 σ 2 1 ∼ N (0, ψ(θ))
. ψ is a positive semidefinite 2 × 2 symmetric matrix (covariance matrix), parameterized by a variance vector θ. We can again regroup the parameters categorized by the group to obtain:

y it = β 00 + β 10 x it fixed effects + b 0i + b 1i x it random effects + it .
(3.12)

Here, we presented a straightforward model with only one explanatory variable. In practice, several explanatory variables can be included in the model. It is not necessary to attribute both fixed and random effects to each of them. We can generally distinguish by their notation: the explanatory variables -which have only a fixed effect -will be grouped together in a matrix denoted X (i.e. fixed effects design matrix), and those -that have both fixed and random effects -will be grouped together in a matrix denoted Z (i.e. random effects design matrix) [Pin, 2000].

Equation (3.12) can be then written as: (3.13) which corresponds to the matrix notation: (3.14) where:

y it = [1 x it ] β 00 β 10 + [1 x it ] b 0i b 1i + it ,
y = Xβ + Zυ + ,
• y is a known vector of observations (e.g. crop yield data), with mean E(y) = Xβ;

• β is an unknown vector of fixed effects;

• υ is an unknown vector of random effects meaning that mean E(υ) = 0 and variance-covariance matrix var υ = G;

• is an unknown vector of random errors meaning that mean E( ) = 0 and variance-covariance matrix var = R;

• X and Z are known design matrices relating the observations y to β and υ, respectively.

Model complexity and overfitting

"With four parameters I can fit an elephant and with five I can make him wiggle his trunk."

-John von Neumann

Many types of models can potentially be used for crop modelling. Their parameters can be adjusted, which leads to different complexity levels, even for one model type. It is crucial to choose an adequate level for each specific application (e.g. statistical crop modelling). This section will present a general idea of why it is important to balance the level of model complexity. The following sections will discuss more detailed methods to find the optimal level of complexity.

Bias, variance, and model complexity

Bias-variance dilemma

Figure 3.3a shows a visualisation of bias and variance. Generally, the bias refers to the error between the average model prediction and ground truth, while the variance refers to the error due to sensitivity to noise in the model prediction for the given datasets. An overly simplified model results in substantial errors because the lack of complexity does not allow to correctly represent the diversity of the target: its output varies little, i.e. low variance, but it is often far from the target, i.e. high bias (Fig. 3.3b). On the other hand, a very complex model (e.g. with a large number of parameters) may also commit errors in the estimation of the target diversity: its output is near the target, i.e. low bias, but it varies a lot, i.e. high variance (Fig. 3.3b). It is known as the bias-variance trade-off or bias-variance dilemma.

Low bias

Low variance

High variance

High bias

Low bias High variance

Model complexity

Prediction error

Low

Model complexity

A too simple model shows high errors on training and testing datasets, which is said to be underfitted (Fig. 3.3b). While a complex model shows a low error on the training dataset and a high error on the testing dataset, it suffers from overfitting (Fig. 3.3b). For example, Fig. 3.4a presents a set of a database in which the blue points are the training dataset, and the orange points present the testing dataset. Here, we propose several tests to perform the regression of y with respect to x: a linear regression (the blue line in Fig. 3.4b), a third-order polynomial (the black curve in Fig. 3.4c), and a fifth-order polynomial (the red line in Fig. 3.4d). Our target is to find a function that fits the training dataset well, but performs well too on the testing dataset. It is shown that with a very simple model (low complexity), the blue line in Fig. 3.4b underfits the data (for both training and testing). In contrast, a fifth-order polynomial, i.e. the red line in Fig. 3.4d, shows clear overfitting (or overtraining). In other words, the red line fits the training dataset perfectly; however, it does not generalize well to the testing dataset. In this case, the reasonable fit is the third-order polynomial, shown in Fig. 3.4c. Its model complexity is more adequate considering the datasets and the underlying model.

Mathematical approach

In mathematical expression, let us consider Ŷ = f (X) the prediction of Y . The average, or the expected value, of the squared difference between the predicted and actual value of Y can also be decomposed into different error terms as follows:

E(Y -Ŷ ) 2 = E[f (X) + -f (X)] 2 = E[ f (X) -f (X)] 2 Bias 2 + E f (X) -E[ f (X)] 2 V ariance + σ 2 Irreducible error (3.15)
The irreducible error cannot fundamentally be reduced by any model. This is because Y is also a function of (the random error term) which cannot be predicted using X. On the other hand, the other two errors, i.e. bias and variance, are reducible errors that represent the inaccuracy in the prediction or difference between f (X) and f (X). To improve the data fitting process of the statistical model (i.e. to select the most appropriate statistical model), we often attempt to minimise these two error terms as much as possible: minimizing the bias and variance errors is the learning process that conducts to the determination of the model parameters.

Overfitting

In general, it is possible that a model performs well with the training dataset but does not perform well with the testing dataset (e.g. as shown in Fig. 3.4d). It means that the model cannot predict the output when dealing with unseen data (data not present in the training dataset). In this case, the model is said to overfit. Several reasons can lead to overfitting. For example, when the model is too complex, this complexity is used to model the noise in the training data (Fig. 3.4d). Or in another case, it is easy to overfit a small dataset (as often in crop modelling tasks, for instance) because the model fails to learn a generalizable mapping. Methods to reduce overfitting will be discussed in the following sections (Sects. 3.5 to 3.7).

Model performance evaluation

Once the adequate model is defined, we need a tool to confirm that the model outputs are acceptable in terms of real data processing, in application mode. In statistics, this step, called model evaluation and validation, is crucial and often done by assessing the model generalization ability since it guides the choice of training method or model. It gives us a measure of the quality of the obtained chosen model. Although very important, this essential step is not always sufficiently considered in the literature, often due to a lack of knowledge in statistical learning.

To evaluate the performance of a statistical learning method on a given data set, we need criteria or metrics to measure how well its predictions match the observation data (i.e. to measure the model quality). Depending on the specific problem, i.e. regression or classification, we will need to use different metrics to evaluate the model. Some of the most commonly used criteria will be discussed in the following for regression and classification problems.

Regression problems

Let's consider n samp samples of yield anomalies y 1 , y 2 , • • • , y nsamp (collectively known as y i or as a vector Y = [y 1 , y 2 , • • • , y nsamp ] T ), each associated with its prediction ŷ1 , ŷ2 , • • • , ŷnsamp (known as ŷi as a vector Ŷ ). The mean of observed data is defined as y = 1 nsamp nsamp i=1 y i .

-The coefficient of determination (R 2 or r 2 or R-squared): is the percentage of the response variation that is explained by a linear model:

R 2 = nsamp i=1 (ŷ i -y) 2 nsamp i=1 (y i -y) 2 = SS Reg SS T ot (3.16)
where SS Reg is the regression sum of squares or also called the explained sum of squares; SS T ol is the total sum of squares (proportional to the variance of the observed data). R 2 is a criterion widely used to compare different regression models. R 2 = 1 means that the predicted (or modelled) values perfectly fit the observed values. A baseline model, which always predicts y, will have R 2 = 0. Models that have worse predictions show a low R 2 . This R 2 value is often considered as explained variances. For instance, R 2 = 0.49 implies that 49 % of the variability of the response variable can be explained by the explanatory variable(s) of the model.

In practice, the use of R 2 should be done with caution since this approach has some limitations. First, R 2 cannot determine whether the coefficient estimates and predictions are biased, which is why the residual plots are often required in addition to R 2 . Also, R 2 is quite sensitive to the number of samples and number of inputs. For instance, R 2 tends to increase with the introduction of new variables, even if they are slightly correlated with the target variable. Consequently, a model with more inputs may appear to have a better fit simply because it has more inputs. In addition, if a model is too complex (e.g. too many predictors or too complex model type), which causes overfitting (Sect. 3.2.2), it produces a misleadingly high R 2 value, and thus it is unreliable in this case. In many cases, the adjusted R 2 ( R 2 or R 2 adj ) is often considered instead of the normal R 2 to deal with the increase of R 2 when extra explanatory variables are added. Among many available adjusting methods, the most used approach is the one proposed by Mordecai Ezekiel [START_REF] Raju | Methodology review: Estimation of population validity and cross-validity, and the use of equal weights in prediction[END_REF]:

R 2 = 1 -(1 -R 2 ) n samp -1 n samp -p -1 (3.17)
where n samp is the number of samples (as above), p is the total number of explanatory variables (not including the constant term).

As shown in Eq. (3.17), the R 2 has been adjusted for the number of predictors in the model. This value increases only if the newly added variables improve the model (i.e. R 2 increases) more than would be expected by chance. The value of R 2 is always less than or equal to that of R 2 .

-The correlation coefficient (COR): is used to measure how strong a statistical relationship is between two variables. Several types of correlation coefficients exist (e.g. Pearson, intra-class, range); we focus here on Pearson's correlation (also called Pearson's R), denoted as COR, a correlation coefficient commonly used in linear regression. Pearson's correlation coefficient is the covariance of the two variables divided by the product of their standard deviations:

COR = cov(Y , Ŷ ) σ Y σ Ŷ (3.18)
where cov is the covariance, σ Y is the standard deviation of Y , and σ Ŷ is the standard deviation of Ŷ . Equation (3.18) can also be expressed in terms of mean and expectation, as follows:

COR = E[(Y -µ Y )( Ŷ -µ Ŷ )] σ Y σ Ŷ (3.19)
where µ Y is the mean of Y , µ Ŷ is the mean of Ŷ , and E is the expectation. The sign of COR indicates the direction of the linear relationship between Y and Ŷ . The absolute value of COR is between 0 and 1 (i.e. ranging from -1 to 1). This value gives us the relationship strength: the larger the number is, the stronger the relationship is. For instance, if |COR| is close to 1, it indicates a strong relationship between Y and Ŷ , while a |COR| close to 0 means this relationship is weak. For a linear model, COR is exactly the square root of the coefficient of determination R 2 (COR = √ R 2 ). Therefore, only one criterion (among COR and R 2 ) should be considered when using linear models.

-The mean square error (MSE): measures the deviation from predictions to observations. This value is calculated as:

M SE = 1 n samp nsamp i=1 (y(i) -ŷi ) 2 (3.20)
In practice, we can also consider the square root of MSE, known as the root mean square error (RM SE = √ MSE), that has the same units as the quantity being estimated. Another common criteria, the mean absolute percentage error

(M AP E = 1 nsamp nsamp i=1 y(i) -ŷi y(i)
), can also be used to measure the accuracy of a forecast system (as a percentage).

-The Akaike information criterion (AIC): is an estimator of prediction error and thereby the relative quality of statistical models for a given set of data. The AIC value of the model can be calculated as follows [START_REF] Akaike | A new look at the statistical model identification[END_REF]:

AIC = 2k -2ln( L) (3.21)
where k is the number of independently adjusted parameters (i.e. estimated parameters), and L is the maximum value of the likelihood function for the model.

In practice, we often compare the AIC values of each model in a set of candidate models to select the adequate model with the smaller AIC value.

In addition to the AIC criterion, one could also consider the Bayesian information criterion (BIC) [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF]. The BIC penalized more heavily more complex models and it is formally defined as: BIC = kln(n samp ) -2ln( L).

Classification problems

Let us take a simple example of a classification problem: classifying between coffee (suitable) and background (not suitable) samples. When running a suitability classification model, our resulting outcome is a binary 0 (not suitable) or 1 (suitable). Then, our resulting classification outcomes can be compared to the actual values of the given observation to judge the performance of the classification model.

Several metrics can be used to evaluate this model, e.g. a confusion matrix (often used to define the measurement metrics), or the AUC-ROC curve [START_REF] Bradley | The use of the area under the roc curve in the evaluation of machine learning algorithms[END_REF][START_REF] Fawcett | An introduction to roc analysis[END_REF][START_REF] Sammut | Encyclopedia of Machine Learning[END_REF][START_REF] Tharwat | Classification assessment methods[END_REF].

-The confusion matrix: visualizes and distinguishes between errors. The general form of a confusion matrix is presented in Fig. 3.5a. Four potential outcomes are used to define the measurement metrics of the classifier [START_REF] Fawcett | An introduction to roc analysis[END_REF][START_REF] Sammut | Encyclopedia of Machine Learning[END_REF]]:

• true positive (TP) represents the number of correctly classified coffee points or samples that are suitable for coffee;

• false positive (FP) describes the number of samples misclassified as suitable for coffee, but actually being background samples, i.e. unsuitable (FP is also known as a type I error);

• false negative (FN) represents the number of samples misclassified as unsuitable for coffee, but actually being coffee samples, i.e. suitable (FN is also known as a type II error);

• true negative (TN) indicates the number of correctly classified background points that are unsuitable for coffee.

The outcomes of the confusion matrix can be used to calculate performance metrics for classification models. The most common metrics are accuracy, precision, recall (or sensitivity), and F1 score (presented in Fig. 3

.5b).

Accuracy is computed as the ratio of correctly classified samples (TP+TN) to the total number of samples (TP+FP+FN+TN): Accuracy = T P +T N T P +F P +F N +T N . Accuracy is commonly used to evaluate model performance. However, one should be careful when using this metric, especially in the case of unbalanced datasets. For instance, if 90 % of cases are true and only 10 % are false. There is a high possibility of obtaining a high accuracy score, e.g. around 90 %. This score does not necessarily imply that the model is good, but it is actually biased toward the dominant class.

Precision is the measure of correctly identified coffee samples over the number of all correctly identified samples: P recision = T P T P +F P . This metric is a good measure to use when we want to focus on limiting false positives.

Recall (or sensitivity) presents the measure of correctly identified coffee samples over the count of actual coffee samples: Recall = T P T P +F N . Using this metric, we can assess how well our model is able to identify the actual true class. It is also useful when we want to focus on limiting the false negatives.

F1 score (or F measure) is the harmonic mean between precision and recall, in which F 1Score = 2(P recision×Recall) P recision+Recall . The harmonic mean is used here instead of the normal mean, as it can better punish values that are further apart. Maximizing the F1 score can thus limit both false positives and false negatives.

-The AUC-ROC curve: is also one of the most widely used metrics for evaluation [START_REF] Bradley | The use of the area under the roc curve in the evaluation of machine learning algorithms[END_REF][START_REF] Fawcett | An introduction to roc analysis[END_REF]. The receiver operating characteristics (ROC) curves are very helpful when we need to identify a good classification threshold as it summarizes how well each threshold performed in terms of true positive rate (T P R = T P T P +F N ) and false positive rate (F P R = F P F P +T N ). In detail, an ROC curve is created by plotting the TPR (or recall) against the FPR for a set of different probability thresholds. For instance, each gray dot in the first ROC curve (ROC 1 plotted in Fig. 3.6) presents the TPR and the FPR for a specific classification threshold. The area under curve (AUC), the area under this ROC curve, is often used to measure the quality of a classification model: the larger the area, the better the performance. AUC ranges between 0 and 1. A value of 0.5 reflects a model that is no better than chance, e.g. the case of the diagonal line in Fig. 3.6. On the other hand, higher values (approaching 1) indicate better model performance: Fig. 3.6 shows that model 1 (corresponding to ROC 1 ) performs better than model 2 (corresponding to ROC 2 ). However, a few drawbacks must be considered before using the AUC-ROC curve [START_REF] Lobo | Auc: a misleading measure of the performance of predictive distribution models[END_REF][START_REF] Peterson | Rethinking receiver operating characteristic analysis applications in ecological niche modeling[END_REF]. For instance, AUC can be misleading in predictive distribution models when background samples are drawn from different background extents [START_REF] Lobo | Auc: a misleading measure of the performance of predictive distribution models[END_REF]. In addition to AUC values, one could also consider a modification of the ROC curve [START_REF] Peterson | Rethinking receiver operating characteristic analysis applications in ecological niche modeling[END_REF] or a calibrated AUC measure [START_REF] Hijmans | Cross-validation of species distribution models: removing spatial sorting bias and calibration with a null model[END_REF] as additional and/or alternative metrics.

Summary

There are plenty of available options to evaluate our model. Depending on the category of the problems (regression or classification), we will need to use different metrics to evaluate our models. In addition, it is important to understand the pros and cons of each metric used so that we can measure correctly how good our models are.

The list of diagnostics presented here are very general and can all be estimated in the learning (or training) database. By increasing the model complexity we can literally improve model performance by getting a better R 2 or AUC values, for example. However, this practice can be misleading as we also need to evaluate our model on the new datasets to check the model quality, choose model complexity, and avoid overfiting.

The training, validation, and testing datasets

One of the main challenges in statistical inference is that the model is set up using a samples database, but it must perform well on new -previously unseen -samples. For that purpose, the overall available database B needs to be divided into three datasets: [START_REF] Ripley | Pattern Recognition and Neural Networks[END_REF].

B = B T rain + B V al + B T est
• The training dataset B T rain is used to calibrate the model parameters once the model structures have been chosen.

• The validation dataset B V al is a sample of data held back from the training dataset, which is used to find the best model (i.e. model selection). For instance, it helps tune the model hyperparameters: choose the more appropriate inputs (i.e. feature selection), determine the number of predictors, find the best model type (linear regression, random forest, neural networks), and determine some training choices. This important dataset is sadly often forgotten, which is often the origin of many misleading results6 .

• The testing dataset B T est is held back from the training and the validation datasets to estimate the true model generalization ability.

The process of partitioning B will be called the "folding" process in the following. For example, the folding choice can be chosen using B T rain = 50 %, B V al = 25 %, and B T est = 25 %.

In general, if the database is big, it is possible to use only training and testing datasets: the training dataset is used to fit the parameters, and the testing dataset is used to estimate the model quality but also to choose the best model. In this case, validation and testing datasets will be similar; therefore, using one testing dataset is acceptable. Nevertheless, the situation becomes different with a small sample size: there is a need for the validation dataset. As in this case, the validation and testing datasets might be different due to the limited sample, and thus results would be misleading. Using only a testing dataset -without a validation dataset -risks choosing the model that best suits this particular testing dataset. This represents a special kind of overfitting, which is not on the model calibration but on the model choice. Thus, using only the testing dataset instead of the testing and validation datasets can be misleading. To avoid this difficulty, we will propose in the following (i.e. Sect. 3.7) a technique that includes a dataset to calibrate the model (training) and another one to choose the best model (validation). The truly independent testing dataset is then used to measure the model generalization ability to process unseen data.

Regularization techniques

Previously, we had discussed several metrics to evaluate the model performance. However, good metrics are often not enough; we also need tools to actually optimise the right complexity of the model to avoid overfitting (Sect. 3.2.2). Regularization techniques can often be used to overcome (or reduce) the overfitting problem, but they can also increase the model interpretability [START_REF] Hastie | Ensemble Learning[END_REF][START_REF] Lecerf | Assessing the information in crop model and meteorological indicators to forecast crop yield over Europe[END_REF]. They can be implemented in multiple ways, for example, by modifying the representation of the data, the loss function, or the learning algorithm [START_REF] Aires | Problèmes inverses et réseaux de neurones : application à l'interféromètre haute résolution IASI et à l'analyse de séries temporelles[END_REF][START_REF] Hastie | Ensemble Learning[END_REF]]. In the following, we will present some commonly used regularization techniques.

Variable selection

In crop modelling, we often need as much information as possible in the inputs to the impact model due to the complex weather-to-yield relationships. However, crop modelling is often characterised by datasets with a very limited number of samples [START_REF] Prasad | Crop yield estimation model for Iowa using remote sensing and surface parameters[END_REF][START_REF] Ceglar | Impact of meteorological drivers on regional inter-annual crop yield variability in France[END_REF]. This small sample size often requires a simple model (i.e. not too high complexity level) to avoid overfitting. Therefore, a selection of explanatory variables is necessary to limit the number of inputs. In other words, this method helps to limit the number of model parameters and thus its complexity. It is part of the parsimony principle [START_REF] Hastie | Linear Methods for Regression[END_REF][START_REF] Gori | Chapter 2 -learning principles[END_REF] that privileges simpler models over more complex ones for the same performance.

In this selection process, we want to have the best combination of predictors (i.e. explanatory variables) statistically. However, this procedure often comes at a very high computational cost as it needs to test many combinations of predictors. Moreover, it is impossible to test all possible models as well as the values of the regularization parameters on each possible combination. More practical methods are commonly considered, for instance, forward selection, backward elimination, or stepwise method [START_REF] Hocking | A biometrics invited paper. the analysis and selection of variables in linear regression[END_REF][START_REF] Xu | Comparison of different methods for variable selection[END_REF]. To perform these methods, an iterative process is often used to provide a hierarchy of the explanatory variables, ranging from a few variables (e.g. in [START_REF] Dinh | River Discharge Estimation based on Satellite Water Extent and Topography: An Application over the Amazon[END_REF]) to many or even about a hundred meteorological variables (e.g. in Mathieu and Aires [2018a]).

We first choose quality criteria (for instance, COR, RMSE, or p-value) to quantify the relevance of keeping a variable as the model's input. Then, a forward selection can be performed as follows [START_REF] Chatterjee | Regression analysis by example[END_REF].

• We begin with a model that contains no variables (i.e. null model).

• The first variable included in the model is the one that provides the best value of the considered criterion, e.g. the greatest correlation or covariance with the response.

• The second one is selected among the remaining predictors as the one that, combined with the first selected input, defines a 2-input model with the best value of the quality criterion.

• This process is repeated until a pre-specified stopping rule is reached (e.g. reaching a desired number of variables), until all the considered predictors are included in the model, or when the testing error actually increases.

In theory, this method first identifies the variables that are most related to the response (e.g. yield anomalies), and then it ensures the selection of the most complementary variables. For instance, two variables that are too correlated would not both be selected.

Once the hierarchy of the explanatory variables is decided, the choice of the number of inputs is also important. We will discuss this choice and how to obtain an optimal number of variables in Sect. 3.7.

In addition to forward selection, backward elimination is also commonly used in variable selection. In contrast to the forward selection, this method begins with the full model (i.e. the model with all explanatory variables) and successively eliminates one at a time. The removed predictors are the ones that do not have a significant effect on the response (often based on the p-value interpretation). The process is terminated when all the variables are significant, or all but one variable has been removed. Depending on the specific application, one should carefully consider using this method as it can be time-consuming by taking into account all the possible explanatory variables.

Another method of variable selection is stepwise method, which alternates between forward and backward techniques [START_REF] Hocking | A biometrics invited paper. the analysis and selection of variables in linear regression[END_REF]. This method is similar to the forward selection method. Nevertheless, at each stage, all variables introduced previously in the model are re-examined. Therefore, a variable that entered the model in earlier stages may be eliminated at the later stages due to its correlations with other introduced variables. The calculations made for the addition or subtraction of variables are the same as two previous techniques (forward selection and backward elimination). The use of this stepwise method has been criticised due to several possible drawbacks [START_REF] Copas | Regression, prediction and shrinkage[END_REF][START_REF] Judd | Multiple Regression[END_REF]. For example, it will not necessarily produce the best model if there are redundant predictors. Models identified by the stepwise method have an inflated risk of making type I error [START_REF] Judd | Multiple Regression[END_REF].

All these methods, by limiting their search, not considering all the possible combinations, are not perfect. However, in practice, we just need to choose one.

Loss function modification

As shown in the example of linear regression (Sect. 3.1.1), the optimal model parameters (or coefficients) are obtained by minimizing the loss function or RSS. Generally, we will first estimate the coefficients based on the training dataset. If there is noise in this dataset, the estimated coefficients will not generalize well on the new unseen data. To avoid this problem, regularization techniques that modify the loss function are introduced [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Gruber | Improving Efficiency by Shrinkage: The James-Stein and Ridge Regression Estimators[END_REF]. Two main types of this loss function-based technique are L2 and L1 regularizations.

L2 regularization (or ridge regression)

Let us consider the same linear regression problem as in Sect. 3.1.1. In L2 regularization, the loss function (Eq. (3.4) in Sect. 3.1.1) is altered by adding a penalty term which is the weighted L2 norm of the weights (β). The optimisation function becomes:

min β ||Y -Xβ|| 2 + λ||β|| 2 , with λ ≥ 0. (3.22)
Here, the penalty term regularizes the model coefficients. Thus, L2 regularization reduces the magnitudes of the coefficients that help to decrease the model complexity. The turning parameter λ is used to control the compromise between overfitting and underfitting: if λ is chosen low (or high), the model will tend to overfit (or underfit). This technique is often used when the explanatory variables are highly correlated (multicollinearity problem). However, L2 regularization has some limitations in variable selection and model interpretability. Because this technique shrinks the coefficients for least important predictors close to zero, but it will never make them exactly zero. Thus, the final model will include all the considered predictors, and it is not good for variable selection (or feature selection), although complementary strategies can be put in place.

L1 regularization (or lasso regression)

Similar to L2 regularization, the L1 regularization uses the L1 norm (absolute values) of the weights in the loss function instead of the L2 norm. The optimisation function can be written as:

min β ||Y -Xβ|| 2 + λ||β||, with λ ≥ 0.
(3.23)

In this technique, the penalty term can better force some of the coefficient estimates to be equal to zero. Therefore, L1 regularization can remove some of the variables and thus can be useful in variable selection too. However, if there are two or more highly correlated variables, L1 regularization will select one of them randomly. It is thus not good for the interpretation of the model. In summary, both L1 and L2 regularization techniques have their own advantages and limitations. The choice of one or another depends on the specific problem. Another alternative to these two techniques is elastic net, which is a linear combination of L1 and L2 regularization [START_REF] Hastie | High-Dimensional Problems: p N[END_REF].

Early stopping

Early stopping techniques [START_REF] Ripley | Pattern Recognition and Neural Networks[END_REF][START_REF] Goodfellow | Deep Learning[END_REF]] also make it possible to limit overfitting when the training is done iteratively, as in neural networks: a higher number of iterations clearly lead to lower training errors, but the model is more sensitive to overfitting. Early stopping aims at stopping the training as soon as the model performance stops improving on the validation dataset. It is a very easy and convenient solution that is always used in the training of NN, for instance.

Dimensionality reduction

Together with regularization techniques, one could also use some dimensionality reduction techniques in variable selection and/or variable reduction. From that, we can reduce the number of variables and thus reduce the model complexity to avoid overfitting. Here, we will introduce one commonly used reduction technique: principal component analysis (PCA) [START_REF]On lines and planes of closest fit to systems of points in space[END_REF][START_REF] Jolliffe | Principal component analysis: a review and recent developments[END_REF]. In general, PCA is a statistical procedure that converts a set of possibly correlated variables into a set of linearly uncorrelated variables called principal components. In other words, it helps us to find a reduced number of variables that will represent our original database in a compressed way, capturing up to a certain portion of its variance depending on the number of new "components" we end up selecting.

In practice, a PCA analysis can be performed in several steps. Let X be a matrix containing the initial data, X ∈ R nsamp×n input .

• First, the initial variable stored in X are standardized such that each has zero mean and unit standard deviation. This step helps to ensure that each variable contributes equally to the analysis. It is optional and can be not good in some cases. Generally, we apply this step when the variables in X have different units.

• The next step is to compute the covariance matrix of the variables in the dataset, and thus compute the eigenvalues and eigenvectors of this covariance matrix.

• Eigenvalues are then sorted in a decreasing order from high eigenvalues to lower, representing a decreasing variance in the data.

• Finally, the transformation (or projection) of the initial (normalised) data onto the reduced PCA space is obtained by multiplying the initial normalised data by the leading eigenvectors of the covariance matrix (also known as the principal components).

In this thesis, the PCA can be used as a regularization technique. In detail, we can feed the statistical models with the dimensional reduced data from PCA. This strategy can help avoiding the multicollinearity problem in linear regression, for instance.

Cross-validation (when having limited data)

"Keep it simple."

Together with regularization and/or dimensionality reduction techniques, one also needs a validation technique to assess how the results of a statistical analysis will generalize to an independent dataset. A simplest method is to divide the overall database into different datasets (shown in Sec. 3.4). However, with a limited number of samples (as often in crop modelling), this division becomes more challenging. In this case, the training process may need every possible data point to determine model parameters [START_REF] Kuhn | Applied predictive modeling[END_REF]. It is thus impossible to keep a significant percentage of the database for the validation and testing datasets. To choose an appropriate model and avoid overfitting, a robust way to measure the generalization ability is necessary, using as few samples as possible.

Cross-validation [START_REF] Allen | The relationship between variable selection and data agumentation and a method for prediction[END_REF][START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF] was introduced as an effective method for both model selection and model assessment when having a small number of samples. In general, cross-validation is a resampling method that uses different portions of the data to test and train a model on different iterations. It can be classified in several types, for example, exhaustive cross-validation (leave-pout cross-validation, leave-one-out cross-validation or simply LOO), non-exhaustive cross validation (k-fold cross-validation, holdout method), or nested cross-validation.

In the following, we will present one commonly used cross-validation method, the LOO, and propose one particular implementation of the nested cross-validation, the leave-two-out nested cross-validation or simply LTO.

Traditional leave-one-out

The LOO method is one common type of cross-validation in which the model uses only two datasets: one to train and another to choose the model and evaluate the result. The main idea of LOO is that given n samp available samples in B, the model is calibrated n samp times (or n samp folds) using (n samp -1) samples in the training dataset B T rain (leaving one sample out). The resulting model is then tested on the left samples (B T est ). There are n samp testing score estimations, one for each sample. In this case, B = B T rain + B T est and B V al is empty. The averaging of these n samp testing scores is expected to be a robust assessment of the model ability to generalize on new samples. However, since no validation dataset is used to select the best model, the choice of the best model may be biased towards this testing dataset [START_REF] Cawley | On over-fitting in model selection and subsequent selection bias in performance evaluation[END_REF]. The chosen model is not independent of the testing dataset, and thus, the obtained testing score is not reliable.

LOO is very useful in many cases [START_REF] Kogan | Winter wheat yield forecasting in Ukraine based on Earth observation, meteorologicaldata and biophysical models[END_REF][START_REF] Li | Toward building a transparent statistical model for improving crop yield prediction: Modeling rainfed corn in the u.s[END_REF] but as described in Sect. 3.4, it is preferable to divide the database into three partitions rather than only two as done under LOO.

Proposed leave-two-out

In the literature, another more complex approach has been introduced to improve the LOO: the nested cross-validation [START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF], also known as double cross-validation or k × l-fold cross-validation, is able to use three datasets: training, validation, and testing. In detail, this approach considers one innerloop cross-validation nested in an outer cross-validation loop. The inner loop is to select the best model (validation dataset), while the outer loop is to estimate its generalization score (testing dataset). In the following sections, we describe one particular implementation of this nested cross-validation (or k × l-fold cross-validation when l = k -1) that we called the leave-two-out nested crossvalidation (LTO), which can improve model selection when the number of samples is low.

Folding scheme

For LTO, we will divide the database into three datasets: training, validation, and testing. Each time we split or partition the dataset, this is referred to as a "fold." Each fold divides the database B into a training dataset B T rain of (n samp -2) samples and a validation B V al and a testing B T est dataset with one sample each. Two samples are considered out of the training dataset instead of one in the LOO procedure. This folding process is presented in Fig. 3.8, with the number of folds n f old = n samp × (n samp -1). This is why this approach is also called k × l-fold cross-validation when l = k -1.

Validation and testing scores

Figure 3.9 illustrates how the LTO evaluation procedure is conducted. In panel (a1), the number of candidate models n mod (represented on the horizontal axis) is defined with a fixed complexity λ of the model. For instance, for the LIN3 model (i.e. LIN model with three inputs) with 12 potential predictors, we obtain n mod = C 3 12 = 220 models. These models are used to perform the yield anomaly estimations. In the vertical axis, for each of the n samp choices of the testing value id test ∈ {1, 2, • • • , n samp }, there are (n samp -1) possible validation datasets and 58 Each choice of the testing value (each id test ) corresponds to a selected best model bm i and two distributions (i.e. probability density functions (PDFs)) for (n samp -1) validation errors and (n samp -1) testing errors, shown in Fig. 3.9a2. These two distributions result in a validation score (blue dot) and a testing score (red dot). The shape of these distributions give the average goodness of fit score and its variance.

Finally, the n samp testing choices give n samp validation and n samp testing scores that form a validation PDF in a blue line, a testing PDF in a red line, and thus the two scores V λ and T λ in Fig. 3.9b.

A pseudo-code is provided below to facilitate the implementation of the LTO procedure in any language.

Generalization ability versus model selection

The process represented in Fig. 3.9 is used to obtain the validation (V λ ) and testing (T λ ) scores from the LTO approach for a given model complexity, for instance, here λ represents the number of inputs. Each different number of inputs (different value of λ) results in different values of V λ and T λ . The V λ and T λ evolution curves obtained for validation and testing RMSE values of yield anomalies for an increasing number of inputs are presented in Fig. 3.10. For simplicity, only validation and testing scores will be discussed since the training error almost always decreases with the number of inputs. At the beginning, both validation and testing errors are large, which suggests that the model is underfitted. Then, it reaches an optimal state, at n input = λ for example. After that, the validation error is smaller but the testing error is bigger when increasing the number of inputs; this is typical for overfitting (Sect. 3.2.2).

Chapter 4

Robusta coffee yield estimation In this chapter, we aim to develop data-driven approaches to investigate the impact of weather on coffee yield. As discussed previously (chapter 1), the main challenge with statistical crop models is that it is based on often limited data. This limitation is no exception to coffee, as exemplified in our case study for coffee-producing states in Brazil as well as districts and provinces in the Central Highlands of Vietnam. To overcome this small sample size issue and avoid overfitting, cross-validation and especially leave-one-out cross-validation (LOO) are widely used in many crop models [START_REF] Kogan | Winter wheat yield forecasting in Ukraine based on Earth observation, meteorologicaldata and biophysical models[END_REF][START_REF] Li | Toward building a transparent statistical model for improving crop yield prediction: Modeling rainfed corn in the u.s[END_REF][START_REF] Laudien | Robustly forecasting maize yields in tanzania based on climatic predictors[END_REF]. In the first approach, we will use a simple multiple regression model with regularization techniques (variable selection and PCA) and LOO validation. We use gridded weather data from the ECMWF re-analysis to identify the most relevant weather information impacting coffee yield. The main objectives are: (1) to propose a simple approach that identifies the most important weather predictors related to coffee yield; (2) to understand the different weather responses of Arabica in Brazil and Robusta in Vietnam; and finally (3) to propose a mathematical framework, focusing on Robusta coffee in Vietnam, to predict the weather-related variability of yield.

The first approach is commonly used in the literature; however, it is often empirical and requires experience. In addition, simple cross-validation can be misleading in several cases. Thus, this chapter will also perform another approach which is a more advanced LTO technique but probably easier to implement. Several models with a varying number of inputs and/or number of potential predictors will be assessed using our proposed LTO and commonly used LOO (without other regularization techniques) techniques. This comparison is performed to investigate the ability in selecting suitable crop models based on the two crossvalidation techniques. Finally, the LTO is exemplified in assessing the Robusta yield estimation model.

In the following, the study area and overview of coffee phenology, particularly the biennial bearing phenomenon, are first presented in Sects. 4.1.1 and 4.1.2. Section 4.1.3 shows the data processing of the yield and weather data. Then, Sect. 4.2 presents the impact model using LOO and regularization techniques, which is the main focus of this chapter. For this model, we will first introduce the general model. We then perform correlation analyses of the yield to different weather variables in different study areas. Results and different aspects of the models are also discussed. After that, a second analysis using a more sophisticated technique (i.e. LTO) is applied for an example of Robusta coffee (Sect. 4.3). More applications using the LTO method will be presented in the next chapter of this thesis (chapter 5). Finally, some conclusions are discussed in Sect.

Materials

Study area

While coffee yield time series were gathered globally, so far, only Brazil and Vietnam have a time record long and suitable for relevant statistical analysis (Sect. 2.1.1). Here, we will focus on five major Brazilian coffee-producing states (shown in Fig. 4.1a) where Arabica is predominant. For Vietnam, the study area covers the high-intensity Robusta coffee production districts (i.e. 20 districts where the coffee planting area is higher than 10 000 ha, shown in Fig. 4.1b) in the four provinces of the Central Highlands.

Overview of coffee phenology

Phenological processes and relationship with weather

As we aim to study the weather impact on coffee yield, it is necessary to understand the phenological processes and how they are affected by the weather. Generally, a coffee plant takes roughly three years to develop, from seed germination to flowering and fruit production [START_REF] Wintgens | The Coffee Plant, chapter 1[END_REF][START_REF] Arcila-Pulgarín | Application of the extended BBCH scale for the description of the growth stages of coffee (Coffea spp.)[END_REF]. It is productive for about 30 years [START_REF] Wintgens | The Coffee Plant, chapter 1[END_REF] but can be more than 50 years, or even up to 80 years with good management. The phenological cycle of coffee trees varies depending on various factors such as the species, variety, weather, and agriculture practices. • Flower bud formation and growth: During this period, the development of a serial bud into a flower bud is largely controlled by plant hormones activated by the change in day length and/or a drop in temperature [START_REF] Descroix | Environmental Factors Suitable for Coffee Cultivation[END_REF][START_REF] Majerowicz | Induction and differentiation of reproductive buds in Coffea arabica L[END_REF]. At the end of the growing phase, these buds enter into a dormant phase until flowering is stimulated.

• Flowering: The dormancy of the flower buds is usually broken by sudden relief of water stress (e.g. the first rains of the wet season, irrigation) and/or a drastic fall in temperature. Water is necessary for this process, but heavy rain could cause an acute form of floral atrophy. The flowers are often pollinated within 24-48 hours.

• Fruit development, including pinhead (a very young fruit) grain expansion, and grain (bean) formation: After pollination, the very young fruits (e.g. pinheads) remain dormant with low water requirements. From the second to third month after blossoming, the fruits start to swell, which increases the water requirements. Thus, this process usually happens during the onset of the rainy season.

• Maturation (fruit ripening): Between the sixth and the eighth months after fertilisation, the fruit reaches maturity, which occurs during the beginning and/or the middle of the dry season since adequate temperature and sunlight promote the ripening of coffee fruits.

The time required from flowering until fruit ripening varies according to various factors (e.g. variety, climate condition, management). In general, it is from 6-9 months for Arabica and about 9-11 months for Robusta.

The biennial pattern

The "biennial pattern" or the "biennial bearing phenomenon" means that the coffee trees produce high yields in alternate years. This phenomenon is often pronounced in unshaded, mature Arabica coffee trees [Wintgens, 2004a;[START_REF] Valadares | Yield and production bienniality of dense coffee plantations under different levels of N and K[END_REF][START_REF] Aparecido | Forecasting of the annual yield of Arabic coffee using water deficiency[END_REF] and particularly pronounced in technified systems in Brazil due to their particular pruning systems. For example, the phenology of Arabica in Fig. 4.2a describes this phenological cycle well, i.e. the entire yield determining factors span over two years: the growing season (vegetative growth) in the first year and the reproduction in the second year. The biennial bearing is a result of competition between vegetative and reproductive processes.

The yield biennial phenomenon is also present in Robusta coffee but less intense due to its mitigation with high input use (e.g. fertilisers, irrigation, and pruning) and alternation of the plagiotropic branches in production [START_REF] Wagner | Crop yield bienniality in groups of genotypes of conilon coffee[END_REF]. It is not related to climatic conditions but is more pronounced where management is deficient [Wintgens, 2004a]. In Vietnam, due to the improvement of management, the yield has steadily increased during the 2008-2017 period. The phenological cycle of Robusta coffee in the Central Highlands (Vietnam) is shown in Fig. 4.2b. In short, in one year, the productive stage and vegetative growth (for the next harvest) overlap in a mature Arabica or Robusta tree.

Therefore, the yield in one year can be impacted by the weather not only in that particular year of harvest but also in the previous year. This effect is particularly pronounced when the coffee trees are exposed to drought or other extremes. Thus, weather effects were analysed for about two years (e.g. 17 to 21 months) prior to the harvest time H 0 . In years without such evident extreme weather conditions, the biennial effect on an aggregated scale such as district or provincial level might be less obvious. However, synchronized phasing of coffee trees' biennial pattern in a region has been documented [START_REF] Valeriano | Estimation of coffee yield from gridded weather data[END_REF] for small locations (e.g. Brazilian municipalities). We have, therefore, also taken it into account as districts and provinces in Vietnam are of similar size.

Data processing

Yield data

As presented in Sect. 2.1.1, the state-level coffee statistics were obtained from the Brazilian Institute of Geography and Statistics (IBGE) from 1991 to 2018 (shown in Fig. 4.3a1). The provincial and district levels coffee statistics were obtained from the General Statistics Office of Vietnam (GSOV) from 2000 to 2018 (e.g. district statistics shown in Fig 4 .3b1). Several preprocessing steps are needed to be done before using these data in the impact models.

Long-term trend definition

To analyse the yield data, the long-term trend was first determined. This longterm trend in coffee yield often describes the changes in agricultural practices (e.g. agricultural management and improved genetics) or indirect factors affecting these practices, such as changes in coffee price [START_REF] Mathieu | Statistical weather-impact models: An application of neural networks and mixed effects for corn production over the United States[END_REF][START_REF] Miao | Responsiveness of crop yield and acreage to prices and climate[END_REF]. Consequently, suppressing this trend from the yield time series allows for removing the impact of non-weather-related factors. A moving average window and a linear function were tested.

Anomaly calculation

Analysis of weather impact on coffee yield was based on the absolute yield anomalies with respect to the previous long-term trend. The coffee yield for year t is defined as y(t), and the long-term trend value is represented by y(t). Both y(t) and y(t) are in 10 3 kilograms per hectare (10 3 kg•ha -1 ). The absolute yieldanomaly a(t) is defined as the percentage variation around the long-term trend:

a(t) = y(t) -y(t) (4.1)
If a(t) > 0, then the yield in year t is higher than average, and vice versa.

Trend identification

The coffee yield time series, shown in Figs. 4.3a1 and b1, suggested that Brazil and Vietnam coffee has different long-term trends. Different trend identification approaches will be applied to the two countries. For Brazil, a piece-wise 5-year moving average window will be used to define the Arabica coffee trend (dotted lines in Fig. 4.3a1) in the two periods : 1991-2002 and 2002-2018 separately. Overall, from 2005 to 2018, the yield anomalies are more stable than in the previous period (see Fig. 4.3a2). The yield data in Vietnam are available for a shorter period compared to Brazil. It is observed that the coffee yield in Vietnam steadily increased from 2000 to 2018 (Fig 4.3b1). A simple linear function was used to define the Robusta coffee trend: y(t) = y 0 + α • t, where y 0 is the yield in the year 2000, and α is the annual rate of increment. Similarly to Brazil's case, the Robusta coffee yield anomalies from 2000 to 2011 show more variability, while fluctuations are lower in the later period, i.e. 2012-2018 (shown in Fig 4.3b2 (shown in Fig 4.3b2). This might be due to the coffee price drop at the beginning of the 21 st century, leading to fewer investments in agronomic management and therefore being more sensitive to weather fluctuations [START_REF] Miao | Responsiveness of crop yield and acreage to prices and climate[END_REF]. In contrast, in the second period, Robusta coffee yield became more stable and weather-independent, most likely due to improving agricultural practices.

Weather data

Monthly mean data

The monthly means of different weather variables were collected from the ERA5-Land and projected into the resolution of the selected administrative levels (e.g. state level in Brazil, provincial and district levels in Vietnam), as introduced in Sect. 2.2.1. Another alternative could be analysed by, for instance, agroclimatic zones. However, in this case, there is a need to gather coffee yield by corresponding zones. This step can mix the coffee information or require yield data at a higher resolution. As a result, we will work at the same resolution levels provided by the coffee yield information.

Direct weather variables

In the literature, statistical approaches often consider precipitation, air temperature (minimum, mean, and maximum), solar radiation, and/or other variables derived from these basic variables [de Oliveira Aparecido et al., 2017;[START_REF] Valeriano | Estimation of coffee yield from gridded weather data[END_REF][START_REF] Kath | Not so robust: Robusta coffee production is highly sensitive to temperature[END_REF][START_REF] Kath | Temperature and rainfall impacts on robusta coffee bean characteristics[END_REF]. This study selects monthly total precipitation (P) and average temperature (T) as explanatory variables for coffee yield. We also considered adding other explanatory variables, e.g. maximum and minimum temperature, solar radiation, and water deficit (at a monthly scale). However, we chose not to include these variables for several reasons: (1) these variables show relatively low correlations to the crop yield anomalies; (2) they are highly correlated to P and T variables, especially in the case of Robusta coffee; (3) it is more reasonable to consider a limited number of explanatory variables when having a limited number of samples to avoid overfitting (Sect. 3.2). The weather data will be evaluated for n months before the harvest season's peak (H 0 ). The value n varies from 1 to 24, depending on the coffee phenology in the study locations. To simplify the application, we evaluated n = 19 for both Arabica and Robusta, Weather anomalies could be considered for coffee yield data (as in Sect. 4.1.3.1). However, the weather trend of the 10 to 20 years is relatively low compared to the inter-annual variations [START_REF] Mathieu | Statistical weather-impact models: An application of neural networks and mixed effects for corn production over the United States[END_REF]. Thus, the long-term trend can be neglected, and the relative anomalies will be estimated based on the longterm average. This average value is computed for each of the n months before the harvest time. In addition, we applied a 3-month moving average centred on the particular month (instead of the monthly data) to reduce the variability at the monthly scale. This variability would introduce instabilities in our analysis due to the short database time length. This approach is also considered a regularization technique (see more in Sect. 3.5).

Bioclimatic variables

We also consider the 19 bioclimatic variables described in O'Donnell and Ignizio [2012] (Tab. 4.1). These variables are derived from the monthly temperature and precipitation values in order to describe biologically meaningful variables such as annual trends (e.g. mean annual temperature, annual precipitation), seasonality (e.g. annual range in temperature and precipitation), and extreme or limiting environmental factors (e.g. the temperature of the coldest and warmest month, and precipitation of the wet and dry quarters). The bioclimatic variables are referred to as Bio for the first year and Bio * for the second year, right before the peak of harvest time (H 0 ).

Weather-to-yield impact model using LOO and regularization techniques

Weather-to-yield impact model

Many regression models could be considered to measure the impact of weather anomalies on coffee yield. As regression models are purely data-driven, a learning database is required to empirically describe the relationships between the inputs and the desired coffee yield. Depending on the available data, models can be chosen from domains such as traditional statistics, Bayesian statistics, machine learning, or artificial intelligence. Yield modelling has been based on linear regression [START_REF] Kouadio | Artificial intelligence approach for the prediction of Robusta coffee yield using soil fertility properties[END_REF], random forest [START_REF] Kouadio | Artificial intelligence approach for the prediction of Robusta coffee yield using soil fertility properties[END_REF], neural networks [START_REF] Mathieu | Using Neural Network Classifier Approach for Statistically Forecasting Extreme Corn Yield Losses in Eastern United States[END_REF], or mixed-effects models [START_REF] Mathieu | Statistical weather-impact models: An application of neural networks and mixed effects for corn production over the United States[END_REF]. Details of these models are presented in Sect. 3.1. Ideally, the complexity of the model is adjusted to the complexity of the relationship under consideration, i.e. the effect of weather on coffee phenology. Thus, a complex regression model would be necessary to model all the processes involved in plant growth and development (described in Sect. 4.1.2). However, in practice, we are often limited to more simple statistical models due to the limited number of samples (see Sect. 3.2). As a result, a simple multiple linear regression method (as described in Sect. 3.1.1) is used here to model the relationship between the observations of coffee yield anomalies a(t) and the weather input anomalies

x i : a = β 0 + β 1 • x 1 + β 2 • x 2 + • • • + β m • x m
where m is the number of input variables.

Figure 4.4 presents the datasets and methods used. In detail, the yield anomaly is induced from the yield for the year t, i.e. yield at the peak of the harvest season H 0 . We then select all potential predictors, called here a D pot set, from the 2-year time series before H 0 , including two sets of variables: (1) n = 19 potential values for each direct weather variable, i.e. 19 realisations for each variable P and T; (2) 38 bioclimatic variables, i.e. 19 variables for year (t -1) (Bio) and 19 variables for year (t -2) (Bio * ). The inputs x i will be chosen from this potential dataset D pot .

Regularization techniques and PCA are then used to avoid overfitting (see more in Sects. 3.5 and 3.6). Firstly, we reduced the number of potential inputs by selecting only four from all potential inputs in D pot , using correlation scores. Secondly, Principal Component Analysis (PCA) was applied over the four selected variables to compress this information into two inputs. By reducing the number of inputs, PCA reduces the number of parameters in the model, therefore regularizing the model structure and reducing the over-training problem. Furthermore, PCA enables the inclusion of more synthetic inputs, with more robust integrative information on robusta coffee. This pre-processing facilitates the work of the regression, so it helps to regularize the problem. We then compared the two approaches: the four selected weather variables and two PCA components to identify the more robust model that could better deal with over-training.

We assessed the Pearson's correlation coefficient (COR) and RMSE (Sect. 3.3) between the estimated yield anomalies (a est ) and the observed yield anomalies (a obs ) as two diagnostics of the model quality. The COR is unitless, and RMSE is in (10 3 kg•ha -1 ).

Finally, the LOO method (Sect. 3.7.1) is used to estimate the model's generalization ability. The same diagnostics COR (i.e. generalization score) and RMSE were considered to measure the model quality.

Correlation analysis

Let us first look at the correlation analysis to understand the coffee yield to weather relation and thus select the four most important predictors.

Yield to direct weather variables

Figure 4.5 presents the correlation between the coffee yield anomalies and the direct weather variables (total precipitation (P) and average temperature (T)) over the two examples: Arabica in Minas Gerais (Brazil) and Robusta in Dak Lak (Vietnam). Only two representative cases are shown here as similar behaviours are found for other Brazilian states and Vietnamese provinces. An opposite relation is observed between these two factors. For instance, while temperature shows a negative correlation (meaning it is beneficial to have lower temperatures), precipitation shows a positive one with the coffee-yield anomaly (meaning that higher precipitation is beneficial) and vice versa. This relation is very consistent in both Brazil and Vietnam case studies. In addition, Arabica (in Minas Gerais) shows a much lower correlation score than Robusta (in Dak Lak).

The same link with the direct weather variables is identified during the maturation period: H -1 → H -4 for Arabica and H -1 → H -3 for Robusta in Fig. 4.5. Both Arabica and Robusta trees benefit from higher temperature rather than increased water availability, which is reasonable because temperature helps boost coffee fruits (as mentioned in Sect. 4.1.2.1).

However, during other phenological processes (i.e. bud development, flowering, and fruit development) there were differences in correlation patterns between Arabica (in Minas Gerais) and Robusta (in Dak Lak). This contrast could be explained by the different weather patterns between the two regions, as shown in the mean temperature and precipitation in Fig. 4.2.

During the bud development (H -17 → H -10 ), Arabica shows a positive correlation with temperature (shown in the black dotted line in Fig. 4.5) and a negative correlation with precipitation (shown in the black solid line in Fig. 4.5). At this stage, the coffee wants to be in the cool season, especially from H -14 → H -10 , which corresponds to the coolest periods in Minas Gerais, thus a positive correlation between the yield and temperature. The correlation between precipitation and yield anomaly at H -6 then increases suddenly and becomes positive. This increase can be explained by the change from the dry to the wet season in Minas Gerais. Also, it is coincident with the fruit development stage, which requires available water. After that, this correlation decreases gradually until the fruits start to become ripe (H -4 ).

On the other hand, in the case of Dak Lak, at the beginning of the flowering (H -11 ), water is needed to break the bud dormancy stage. In the Central Highlands, this is commonly done using irrigation before the onset of the rainy season. As a result, the presence of precipitation is not so significant for the coffee yield. At a later period, H -4 → H -9 adequate rainfall during the rainy season becomes more important as it corresponds to the fruit development stage which is the stage of highest water demand. This explains the positive correlation with precipitation during this stage. Additionally, the highest positive correlations were found at the end of the rainy season of the previous year (H -13 ) and the peak of the rainy season of the year of harvest (H -4 ) with correlations of 0.50 and 0.55, respectively. This points out that the length of the rainy season is important for successful bean filling and hence final yield. Negative correlation was found for temperature (-0.45 and -0.40, respectively), pointing out that too high temperature can accelerate fruit growth and ripening to the detriment of yield.

Yield to bioclimatic variables

Similarly, Fig. 4.6 presents the correlation between the coffee yield anomalies and the bioclimatic variables. These scores generally are much lower than the ones of direct variables. Also, there are differences between Arabica and Robusta and between the first and second year before the harvest time (i.e. Bio in dotted lines and Bio * in solid lines) on the correlation between yield anomalies and weather variables. In a previous study, Bunn et al. [2015b] stated that among the bioclimatic variables, Bio10, Bio5, and Bio8 are important for Arabica's geographical distribution; while the Bio2, Bio5, Bio7, and Bio15 are ranked high for Robusta. This is in agreement with our results in Fig. 4.6 showing a high correlation with Bio10 * , Bio5 * , and Bio8 * for Arabica (but not for Bio10, Bio5, and Bio8). It is also found that Bio2 * ; Bio5 * , Bio7 * , and Bio15 * are highly correlated to the Robusta yield anomalies.

As mentioned in Sect. 4.1.3.2 and Fig. 4.4, our approach is more directly related to the yield in a more quantitative way: we consider the two sets of bioclimatic variables Bio and Bio * (corresponding to the two years before the harvest time). On the other hand, Bunn et al. [2015b] studied the link between the climate variables (Bio) and the present coffee production occurrence. As a result, there are differences in the significance of Bio or Bio * with Bio variables in the study of Bunn et al. [2015b].

Variable selection

Since both direct weather and bioclimatic variables are derived from the precipitation and temperature data, there is a link between the two datasets. For example, bioclimatic variables mainly consider precipitation and temperature during the wettest/driest or warmest/coldest month/quarters. These variables are similar to the direct weather data, for instance, at H -3 , H -6 , H -16 , or H -18 for Arabica; H -3 /H -4 and H -8 /H -9 for Robusta.

Similar results between the two sets of variables are found for both Arabica and Robusta cases. Figure 4.6 showed that the most important variables for Arabica are the precipitation of the warmest quarter (Bio18), which is the same as the direct weather data precipitation at H -6 (solid black line in Fig. 4.5). For Robusta, the precipitation of the wettest month (Bio13) shows a high correlation with the yield anomalies. This parameter is similar to the high score of precipitation at H -3 or H -4 .

Each weather dataset has its own advantages. For example, using the direct weather data, we can consider not only precipitation and temperature but also other parameters (if needed) like evaporation and surface net solar radiation. Besides, with the direct weather data, we can study a significant period related to the coffee phenology instead of the effect of the whole year. On the other hand, some bioclimatic variables, such as the mean diurnal range (Bio2), seem to be important for the coffee yield.

In this application, we finally select only direct weather variables as the potential predictors for several reasons. For example, higher correlation scores are found for the direct weather variables, which suggests that more information can be exploited from this variable. In addition, considering the available yield database size, it is more reasonable to consider a limited number of explanatory variables to avoid overfitting [START_REF] Lecerf | Assessing the information in crop model and meteorological indicators to forecast crop yield over Europe[END_REF][START_REF] Dinh | Climate change impacts on robusta coffee production in vietnam[END_REF]]. In the following analysis of this chapter, we will also focus only on Robusta in the Central Highlands of Vietnam because (too) low correlation scores are found for Arabica (as shown in the example of Minas Gerais in Figs. 4.5 and 4.6). These low scores could be due to the coarse resolution (i.e. state level) that the model has worked on.

Focused on 20 districts in the Central Highlands of Vietnam

At the district level, the variables (i.e. direct weather variables) with the highest correlation vary slightly among the six considered districts of the Dak Lak province (Fig. 4.7). Due to the climatic variation among districts, it is necessary to use multiple variables as inputs to the yield estimation model. Here, we selected the four most important weather variables, i.e. the precipitation and temperature in the two most important months. Figure 4.7 shows the two months with the highest correlations between weather variables and yield anomalies. These selected months vary from district to district. Overall, the most important month (i.e. shown in Fig. 4.7) ranges from the sixth to the third month before the harvest (H -6 to H -3 ), corresponding to the end of fruit development stage (the bean formation). The second most important month (i.e. shown in Fig 4 .7) is related to the vegetative growth and bud development stages, i.e. H -17 to H -12 .

Results of the impact model Dak Lak district-level case study

Table 4.2 shows the learning and generalization scores, or the correlations between observed and estimated yield anomalies over the Dak Lak districts, for the two input configurations: the four selected weather variables and the two PCA components. The correlation for the model using four selected variables is much higher over the training dataset than over the generalization (Tab. 4.2). This indicates that over-training is occurring when using too many inputs. For instance, in Buon Ma Thuot, the correlation between a est and a obs reduces approximately by 0.4. On the other hand, the learning and generalization scores are more similar when using PCA, which means that the over-training is reduced. Besides, we obtained better generalization scores when using two PCA components instead of the four selected variables. A more realistic generalization correlation/score of 0.55 is obtained for Robusta coffee in Buon Ma Thuot, and the average correlation of 0.49 for the whole province of Dak Lak. It confirms that PCA is a good regularization technique, which will be used in the subsequent models. The generalization correlation ranges from about 0.4 to 0.6, which shows considerable differences among districts. Different environmental conditions can explain these differences, not only between but also within the provinces. observed and estimated yield anomalies for generalization over the selected districts in the Central Highlands. The explained variance (in %), which equals the square of the correlation coefficient, is also indicated (see in Fig. 4.8b). Figure 4.8 shows that the model works best over the Dak Lak districts, followed by Lam Dong, while lower generalization scores are found for Kon Tum and Dak Nong districts. Weather explains 25 % of the variation in Robusta coffee yield anomalies in Dak Lak on average and about 16 % for the whole Central Highlands. These values could appear to be low, but it should be kept in mind that the objective here is not to forecast the whole coffee yield anomalies but only their weathersensitive part. Our experiments show that between 16 and 25 % of the anomalies can be explained by weather, which is very significant for the farmer revenues.

All selected districts in the Central Highlands

Yield estimation

District level case study in Dak Lak Figure 4.9 presents the observed and predicted coffee yield anomalies, as well as the corresponding coffee yield (in kg•ha -1 ) time series and its uncertainties in one Dak Lak district (Cu M'gar). As mentioned above, we obtain a correlation of 0.6 between the observation of yield anomalies a obs and its estimation a est in Cu M'gar. The model can explain 36 % (= 0.6 2 ) of the variation in coffee yield anomalies in this district. The correlation between y obs and y est (0.72) is certainly higher than the yield anomalies correlations due to the use of the trend that is here known without uncertainty (Fig. 4.9). The prediction seems to work better in the 2000-2011 period compared to the 2012-2018 period.

Assessment of the model quality

Besides the yield anomalies, it is also interesting to compare the estimated coffee yield (y est ) and the observed yield or the real data (y obs ). The yield forecasting for year t is estimated as y est (t) = y(t) + a est (t), where y(t) is the yield trend (Sect. 4.1.3.1) and a est (t) is the forecasted yield-anomalies.

In addition, to assess the model quality, we compared the two cases: the average state given by the long-term trend (i.e. no weather information) and our weather impact model estimation. The average state is the most straightforward guess/prediction without weather information, resulting in a typical year with no weather anomaly. We compared this simple guess to our model estimation to see how much we gain by using weather information.

Figure 4.10 presents the results of the yield anomalies (in 10 3 kg•ha -1 ) and production (in 10 6 kg) over the selected districts in 2011 of the three cases: i) the true observation, ii) the average state with no weather information, and iii) our weather impact model. As indicated in panel a2 (Fig. 4.10), the true observation of the total production over the selected districts is 935×10 6 kg in 2011. For this year, our weather impact model shows a similar estimation compared to the true observation (i.e. 926×10 6 kg). Although our model seems to overestimate one Dak Lak district (with an absolute difference of about 16×10 6 kg), the integrated prediction is closer to the true observation.

Sensitivity to the model spatial scale

Each region/district differs in terms of elevation, weather, and perhaps even agricultural practices. Thus, calibrating a localised model could be beneficial if the model can take into account these specificities. However, this is not always possible as sufficient yield data are often unavailable or the weather data is too coarse; therefore, evaluating the differences between scales is of interest. Also, when working at a large scale, we can compensate for all small region uncertainties, thus gaining better results. Hence, working at each spatial scale has advantages and disadvantages, depending on the problems and applications. Therefore, it is necessary to study this scaling aspect and identify which spatial scale is more optimal to work on.

In the case of Robusta in the Central Highlands, the yield data are available at both provincial and district scales. Therefore, we analysed our regression model at both scales to better understand the following questions: (1) at which scale can yield be better forecasted, (2) do we gain more when we go to higher spatial resolution. We first set up the two models at the provincial and district levels. However, to perform a comparison between the models, we applied the provincial model to the district scale to compare it with the district model, and we also integrated the district model prediction at the provincial scale to compare these results to the provincial model.

Figure 4.11 shows the generalization score of different models when applying to different spatial scales (i.e. province and district). Firstly, considering the comparison at the provincial scale: Fig. 4.11a1 is the result of the provincial model, and Fig. 4.11b1 is induced by integrating the district model at the provincial scale. The average generalization score obtained over the Central Highlands is a little higher in Fig. 4.11b1 compared to Fig. 4.11a1, i.e. 0.51 and 0.44, respectively. These scores indicate that the integration of the district model is more interesting than the provincial model.

Similarly, we consider the two cases at the district scale: applying the provincial model to the district scale (Fig. 4.11a2) and using the district model (Fig. 4.11b2). Although some Lam Dong districts show high generalization scores (0.6 or 0.7), the average score over the whole study area is very low (0.25) in Fig. 4.11a2, i.e. the lowest one among four tested cases. Meanwhile, when using the district model, the overall score is much higher (0.39). In short, as expected, for Robusta coffee in the Central Highlands, it is more advantageous to use the model calibrated at the district scale as it can take advantage of the local specificities, even if it is more ambitious.

Discussions of the Robusta analysis

Correlation patterns in relation to coffee phenology

As shown in Sect. 4.2.2.1 for Dak Lak province, the data suggest that Robusta coffee yield (H -2 → H -1 ) benefits from higher temperature during the fruit maturation stage, while less precipitation appears to be required. This might be related to the temperature seasonality, with a high average temperature in March/April and a low average temperature in October/November (Fig. 4.2). As the fruit maturation period coincides with the coolest season of the year, the positive correlation between yield and temperature of that season could indicate that too low temperature slows the maturation rate to the detriment of yield; hence higher temperatures are more beneficial (Sect. 4.1.2.1). Too much water (particularly strong/intense rainfalls), on the other hand, could lead to pre-mature fruit dropping. In the period from H -8 → H -3 , which coincides with the rainy season and the fruit development stage, higher precipitation positively affects yield, while higher precipitation is accompanied by lower temperature. At the beginning of the flowering stage (H -11 ), there is a need for water to break the bud dormancy stage. Coffee farmers in Dak Lak generally start irrigation during this time to induce flowering before the rainy season starts. As a result, the precipitation during that period has not been identified as important for coffee yield.

Different choices of weather inputs for different districts

The districts differ regarding the coffee planted area, mean elevation, and climate (precipitation and temperature). As a result, we identified different constellations of the two most important months where weather affects yield (in Fig. 4.7). Yet they all coincide with the end of the fruit development stage and vegetative growth (bud development stage), highlighting the particular importance of weather during these two phenological stages in the case of Robusta coffee in Vietnam. The correlations indicate that the previous year's prolonged rainy season benefits vegetative growth and hence potential yield, while sufficient rainfall during fruit growth yet low rainfall and higher temperatures during fruit maturation and harvest is more favourable. These moments are of interest to test the yield forecasting, whereby the district coffee yield anomalies could be forecasted from three to six months before the harvest, depending on the properties of the study districts (elevation, climate). In the areas at higher elevations (i.e. 600-1000 m), the weather at six months ahead of the harvest is the most important variable, while at lower elevations (i.e. 200-600 m), the coffee yield anomalies are impacted more by the precipitation and temperature at H -3 . However, the accuracy of such a yield forecast and its value for decision-making would need to be tested in more detail and depend on seasonal climate forecasts for the different districts of interest.

Comparison with previous studies

Our results for Robusta coffee in Dak Lak agree with [START_REF] Kath | Not so robust: Robusta coffee production is highly sensitive to temperature[END_REF] who showed that rainfall has a positive effect while temperature has a negative effect during the growing season (Mar-Sept). One advantage of our approach is that we can point out the month with the highest correlation, while [START_REF] Kath | Not so robust: Robusta coffee production is highly sensitive to temperature[END_REF] only looked at two seasons (flowering and fruit growth). As phenology dynamically responds to prevailing weather conditions, the key moments that affect yield vary in space. For instance, Fig. 4.5 shows a positive correlation score from H -9 to H -3 (i.e. the growing season), with the highest score at H -5 for Robusta coffee in Dak Lak. However, the different districts have different key moments where weather most affects yield. In our example of Dak Lak, we found a low correlation score of precipitation at H -11 /H -10 (i.e. the beginning of the flowering season) (Fig. 4.5), while [START_REF] Kath | Not so robust: Robusta coffee production is highly sensitive to temperature[END_REF] concluded that rainfall has a negative effect on Robusta coffee yields during the flowering season (Jan-Mar). It is unclear whether this difference is due to the different spatial scales of used weather data or due to the underlying yield data set.

The consideration of the biennial pattern

Because Robusta coffee in Vietnam is intensively managed [START_REF] Amarasinghe | Toward sustainable coffee production in Vietnam: More coffee with less water[END_REF][START_REF] Byrareddy | Sustainable production of robusta coffee under a changing climate: A 10-year monitoring of fertilizer management in coffee farms in Vietnam and Indonesia[END_REF], we expected that the biennial pattern might not be relevant to the coffee yield data. However, the weather variables of the previous year before the harvest did affect the coffee yield, as indicated in the correlation analysis (shown in Fig. 4.5). For instance, both precipitation and temperature at H -14 and H -13 are highly correlated to Robusta coffee yield anomalies in Dak Lak. Interestingly, an opposite relation is found between H -13 and H -1 , although both refer to November. As for H -13 , the correlation (between temperature and yield anomaly) is negative, and for H -1 , the correlation is positive, indicating that H -13 negatively affects the yield of the previous year in favour of the yield of the next year. Although it is not clear about the need to consider the biennial pattern in Robusta coffee yet, our study showed that it is important to look at both the current and previous years.

Limitations

This study highlights the importance of selecting spatially dynamic weather variables for optimal modelling. This needs to be atomised to apply this model approach at a larger scale. Therefore, a practical automatic selection procedure needs to be developed. Secondly, given the relatively high weather sensitivity already detectable with simple linear regression using a small yield database, there is a need to consolidate high-quality yield data from all coffee growing regions to enable the development of more complex data-driven approaches and elucidate key weather-sensitive coffee phenological stages and related yield forecasting skill. We expect to obtain better forecasting scores in regions where coffee is rainfed and thus is more sensitive to the weather.

Weather-to-yield impact model using LTO

Previously, we considered a data-driven approach with regularization techniques, dimensionality reduction, and the LOO method to avoid overfitting the regression models. In this section, we will perform the statistical yield modelling of Robusta coffee in Cu M'gar (Dak Lak, Vietnam). The goal is to find a model that makes the most robust predictions of crop yield anomaly as a function of weather variables. We first assess several models (with a varying number of inputs or number of potential predictors) to find the appropriate model choices using both LOO (Sect. 3.7.1) and LTO (Sect. 3.7.2) approaches.

Statistical yield models

The model Similar to Sect. 4.2.1, a multiple linear regression is used here to model the relationship between the yield anomalies and the direct weather variables:

a = β 0 + β 1 • x 1 + β 2 • x 2 + • • • + β m • x m ,
where m is the number of input variables. In this case, m inputs will be selected from a set of direct weather variables: 2 × 19 potential values of P and T (as presented in Sect. 4.1.3.2).

Model selection

Model selection is the process of selecting one model -among many candidate models -that generalizes best [START_REF] Hastie | Model assessment and selection[END_REF]. This process can be applied across models of the same types with varying model hyperparameters or across different model types. This section investigates some practically important factors of the model selection:

• Number of inputs: The inputs are variables that are necessary for model execution through algorithms. The inputs are selected among the potential predictors. We often have a big set of potential predictors (e.g. all-weather variables during the crop growing season), but we select only some variables from this set as the model inputs. The number of inputs defines the model complexity: the higher the number of inputs is, the more complex the model is (supposing that other factors are fixed).

• Number of potential predictors: The potential predictors (i.e. potential explanatory variables) here refer to all possible variables that can potentially impact the yield. Here, we consider 38 weather values for Robusta coffee, but these numbers could be much larger. For instance, in addition to selected weather variables, we could consider other variables (e.g. water deficit, soil moisture), and agroclimatic indices (e.g. degree days, frost-free period [Mathieu and Aires, 2018a]). In addition, monthly variables are used here, but weekly or daily variables could have been considered. Therefore, establishing the list of potential predictors is not fixed in the model selection: it is a crucial modelling step preliminary to any input selection [START_REF] Ambroise | Selection bias in gene extraction on the basis of microarray gene-expression data[END_REF][START_REF] Hastie | Model assessment and selection[END_REF].

Model assessment

Finally, the LOO and LTO methods are used to estimate the model's generalization ability. Here, COR and RMSE are also used to measure the quality of the model.

Yield model selection

We first investigated the model choice by varying the number of inputs. In this example, the number of potential predictors is fixed at 18 (n pre =18). The number of inputs is chosen from one to six, as shown on the horizontal axis in Fig. 4.12. We used the LOO and LTO procedures to compute the corresponding training, validation, and testing RMSE values. The results from the LOO procedure (in Fig. 4.12a) tell us that a model with more inputs is preferable: both training and testing RMSE values decrease with the increase in the number of inputs. In the LTO case, the training and validation RMSE values decrease with the model complexity, similar to the training and testing errors in the LOO procedure. This similarity is because the LTO validation dataset has the same role as the LOO testing dataset: to find the best model! In the case of a too simplistic model, i.e. the LIN model with one input, underfitting occurs as the errors are high in the training, validation, and testing datasets (shown in Fig. 4.12b). These errors decrease gradually with the number of inputs, i.e. from one to three. However, the testing errors do increase when the model has more than three inputs. The LTO procedure indicates that a simple model -with only three inputs -is optimal.

Figure 4.13 shows the RMSE values of the predicted Robusta coffee yield anomalies for the LIN models, with the number of potential predictors ranging from 5 to 38 (on the horizontal axis). These values are computed using the LOO and LTO procedures for the training, validation, and testing datasets. Several models have been tested; we presented here a particular example of the LIN5 model, which is the linear regression model with five inputs. These inputs are selected among the considered potential predictors. For instance, for the LIN5 model with six potential predictors, LOO and LTO aim at choosing five inputs from {P N ov(t-1) , P N ov(t) , T M ar(t) , T Jan(t) , T M ay(t) , and P Oct(t-1) }.

The LOO procedure suggests that the more potential predictors the models have, the better the results are. Both training and testing RMSE values decrease gradually (Fig. 4.13a) with the increase in the number of potential predictors for LIN5 models. On the other hand, the same behaviour is observed for the LTO procedure in Fig. 4.13b; the testing errors show an opposite trend to the training/validation errors and gradually increase with the number of potential predictors. The LTO procedure indicates that a simpler model with fewer potential predictors is more appropriate. This conclusion makes sense since it is inappropriate to use a very complex model (as the LOO model choice) when having a limited sample.

The LOO procedure is actually misleading because it could encourage us to choose a model that overfits the data: the same testing dataset is used to choose the best model and to assess the generalization ability. If the modellers select the best model based on information from the LTO procedure, they are less likely to choose an overfit model. As in this case, they choose the model on the validation dataset and assess its generalization score on an independent testing dataset.

In short, considering the limited information in the available database -that is used to train, select the model, and evaluate its quality -it is not possible to use more than a very simple and limited model. Therefore, for this 19-sample coffee yield modelling case, using a simple LIN model is better than a complex one.

Yield anomaly estimation

The previous section showed that the LTO procedure allows us to choose a reasonable model, simple enough, with fewer inputs and potential predictors. Thus, the crop yield estimations of the LTO method will be assessed here to see how good the selected model (LIN3 model with three predictors) is. The final model includes {P N ov(t-1) , P N ov(t) , T M ar(t) } and these selected variables coincide with the key phenological phases of Robusta coffee. For example, there is the need for a dry period for the buds to develop into dormancy at the end of the development stage, i.e. Nov(t-1) [START_REF] Schuch | Flowering, ethylene production, and ion leakage of coffee in response to water stress and gibberellic acid[END_REF]. Therefore, P N ov(t-1) impacts the buds directly and thus the potential yield. Similarly, the fruit maturation stage (Nov(t)) benefits from weather conditions with less precipitation. At the beginning of the fruit development period (Mar(t)), a too low temperature slows the maturation rate to the detriment of yield, while a higher temperature is beneficial [Wintgens, 2004a].

Figure 4.14 presents the estimated yield anomalies time series for Robusta coffee in Cu M'gar from 2000 to 2018. The estimation (a est in the dashed line) describes quite well the observations (a obs in the solid line) with a correlation of 0.57. With precipitation and temperature variables, the selected model is able to identify many extreme years (e.g. 2005-2009, 2010, 2011) or a decreasing trend from 2011 to 2015. Also, the correlation score means that the model can explain more than 30 % (0.57 2 ) of the variation in coffee yield anomalies, which is in agreement with the result obtained in Fig. 4.9. This value is reasonable as the weather is among several factors (e.g. prices, socio-technical factors, managerial decisions) affecting coffee yield [START_REF] Miao | Responsiveness of crop yield and acreage to prices and climate[END_REF][START_REF] Kc | How climatic and sociotechnical factors influence crop production: a case study of canola production[END_REF][START_REF] Liliane | Factors affecting yield of crops[END_REF]. It is possible to apply the resulting statistical crop yield model to future climate simulations and then study the impact of climate change on coffee [Bunn et al., 2015b;[START_REF] Craparo | Coffea arabica yields decline in tanzania due to climate change: Global implications[END_REF][START_REF] Läderach | Climate change adaptation of coffee production in space and time[END_REF].

Discussions and conclusions of the chapter

This chapter explored how the relationship between phenological processes of the coffee plant (i.e. bud formation and growth, flowering, fruit development, and maturation) and weather affects coffee yield. Each phenological stage requires specific weather characteristics for optimal yield. Due to the perennial nature and biennial pattern of the coffee plant, climatic effects in the previous year can affect the following year's yield. This understanding is essential in studying weather patterns that affect coffee yield. Although analyses were performed for coffee regions in Brazil and Vietnam, we finally focused on Vietnamese Robusta coffee at district and provincial scales because of its high sensitivity to weather.

Two statistical approaches were performed to estimate the weather impact on Robusta coffee yield. The first approach uses a simple regression model with regularization techniques and LOO. The analysis indicates that precipitation and temperature during (1) vegetative growth of the previous year and (2) bean formation are most influential in determining coffee yield anomalies. The exact timing of these variables varies among coffee districts due to the difference in elevation, climate, and/or agricultural practices. Overall, the model can forecast the coffee yield anomalies from three to six months before the harvest, six months for the districts located at a higher elevation with cooler temperatures and, therefore, slower development rates. Nevertheless, the accuracy of such a yield forecast and its value for decision-making would need to be tested in more detail and depend on seasonal climate forecasts for the different districts of interest. Our findings indicated that the weather could explain up to 36 % (from 16 to 25 % on average) of the variation in Robusta coffee yield anomalies in the Central Highlands. This value is in line with the literature [START_REF] Ray | Climate variation explains a third of global crop yield variability[END_REF][START_REF] Craparo | Coffea arabica yields decline in tanzania due to climate change: Global implications[END_REF]. The remaining variability is rather large and may be explained by non-climatic factors (e.g. prices, sociotechnical factors, managerial decisions, or political and social context). It could also come from climate; however, the model would require more detailed variables (e.g. at a daily scale) or more samples to go into greater details of the climate-crop-yield relationship. In addition, this value is already significant as it means that weather can affect up to one-quarter of the farmer's revenue. Some Lam Dong districts can show higher weather sensitivity. Other coffee districts in Gia Lai or Dak Nong are not so sensitive to weather, which could partially be due to the adequacy of these regions to the coffee production and/or higher stability of the weather in these regions.

Identifying an appropriate statistical model for a relatively complex crop like coffee is not easy, particularly as there are limited observations (only about 20 years of coffee yield). This explains why using regularization strategies such as favouring simple models (e.g. linear regression) over complex ones and applying good quality assessment diagnostics are required. In particular, we used the LOO method to avoid potential over-training when a limited number of samples is available. However, this approach is technical and requires experience to master. Although, results from LOO can be misleading in some cases if the analysis was not done probably, without regularization, for instance.

In this chapter, we also perform another statistical approach, which is a particular form of nested cross-validation approach named as LTO method, to calibrate, select, and assess competing models. With the example of Robusta coffee in Cu M'gar, our results show that LOO (when used alone, without other regularization techniques) can be misleading because it uses only one dataset to choose the best model and estimate its generalization skills simultaneously. This is a true problem as LOO is one of the main statistical tools to obtain crop yield models. We, therefore, consider the LTO method that helps select the best model by using two datasets that are independent of the training dataset: first, the validation dataset is used to select the best model form or complexity, and then the test dataset is used to assess the model performance independently. In our case study of Robusta yield prediction, LTO shows that only very simple models can be used when the database is limited in size. It was shown that considering the available historical yield record, the best statistical model can explain about 30 % of the coffee yield anomaly variability, which is comparable with our first approach. The LTO implementation proposed here is very general and can be applied to any statistical crop modelling problem when the number of samples is small and a large number of potential predictors are available. More applications of this LTO In chapter 4, we saw an application of LTO for Robusta yield estimation. This proposed LTO technique is very general and can be applied to any crop and is not limited to coffee. For instance, all crops presented in Sect. 2.1.2 (e.g. French oat, maize, or EU wheat) can be investigated. This chapter focuses on French grain maize, which is very different from coffee for several reasons. Firstly, it is an industrialized agriculture on larger landholdings with France being a leader for over 50 years. Secondly, grain maize is an annual crop and is therefore much simpler in terms of phenology. Because the crop is annual, every year with yield data is independent from previous years and the growing season is only 6 months rather than 12 months.

In the following, the LTO approach will be applied over 96 French departments for several purposes: to obtain a reliable assessment of the model generalization ability, to compare the performances of different predictive models, and to determine the optimal complexity of the statistical grain maize models. In this application, with (a little bit) more yield data and a simple crop, we expect to have models of better quality. There is also the potential to test more complex models like NN. Finally, this application helps to illustrate the flexibility and automation of the LTO technique, which can be potentially applied to any agricultural crop. The following sections of this chapter will (1) introduce the materials and databases used for statistical crop modelling, (2) describe the model, (3) evaluate and select the "best model" by using LOO and LTO approaches, and (4) assess the seasonal yield anomaly forecasts.

This thesis chapter is the subject of a publication in Geoscientific Model Development [START_REF] Dinh | Climate change impacts on robusta coffee production in vietnam[END_REF].

Materials

Overview

Grain maize (Zea mays) is among Europe's most common annual crops. France, our study region, is the largest grain maize producer in Europe [EUROSTAT, 2021]. The study area has improved a lot in agricultural management including irrigation practices since 1960 [START_REF] Siebert | A global data set of the extent of irrigated land from 1900 to 2005[END_REF][START_REF] Schauberger | Yield trends, variability and stagnation analysis of major crops in France over more than a century[END_REF][START_REF] Ceglar | Time-varying impact of climate on maize and wheat yields in France since[END_REF]. Although the sowing time varies for different regions [START_REF] Olesen | Bibliography Changes in time of sowing, flowering and maturity of cereals in europe under climate change[END_REF], the average growing season of French grain maize ranges from April to September [START_REF] Ceglar | Linking crop yield anomalies to large-scale atmospheric circulation in Europe[END_REF]Agri4cast, 2021]. Many previous studies showed that grain maize yield is sensitive to weather conditions [START_REF] Ceglar | Impact of meteorological drivers on regional inter-annual crop yield variability in France[END_REF][START_REF] Ceglar | Linking crop yield anomalies to large-scale atmospheric circulation in Europe[END_REF][START_REF] Lecerf | Assessing the information in crop model and meteorological indicators to forecast crop yield over Europe[END_REF]. Therefore, we will analyse the relationship of maize yield to meteorological variables during the 6-month growing season.

Data processing

Yield data

As presented in Sect. 2.1.2.3, the French grain maize data (area, production, and yield) on the regional level (i.e. department) were collected from the Agreste website for a period of 22 years (from 1989 to 2010). Here, we have modelled the yield of grain maize in 96 French departments (Fig. 5.1). Some specific tests (in Sect. 5.2) will focus on ten departments (as presented in Fig. 5.1d) where the average grain maize production is higher than 4 × 10 5 tons (or the area is higher than 40 000 ha).

The grain maize anomalies are calculated by removing the long-term yield trend, as done in the coffee case study (Sect. 4.1.3.1). A 10-year moving average window is used because the trend is slightly more complex than a linear trend.

Weather data

Similar to the coffee case study of the previous chapter, we will use the direct weather variables as inputs (Sect. 4.1.3.2). Six months of growing period (i.e. from sowing to harvest) will be studied, resulting in 2 × 6 potential predictors: P and T from April to September (P Apr , P M ay , • • • , P Sep and T Apr , T M ay , • • • , T Sep ).

Statistical yield models

The statistical models measure the impact of weather on crop yield anomalies, which can be denoted as: a(t) = f w (X), where f w is the parametric (or nonparametric) statistical model, w stands for the model parameters, and X is the set of weather inputs {X i for i = 1, 2, • • • , n input }. In this chapter, two statistical models are considered: (1) multiple linear regression (LIN) and (2) neural networks (NN). Detailed descriptions of these two models are presented in Sect. 3.1.

Again, two diagnostics are considered, in this case, to measure the quality of the yield anomaly estimations: COR (unitless) and RMSE. The RMSE includes systematic and random errors in the model. The RMSE unit is the same as a(t); RMSE = 40 represents an anomaly error of 40 %.

Model selection and evaluation

As in Sect. 4.3, we investigate different factors relevant to model selection, for example, the number of inputs and the number of potential predictors. In addition, this grain maize application will also include the model types. We perform the selection using two types of models (presented in Sect. 5.1.3) with different levels of complexity. For example, with n input inputs, a simple LIN model usually requires (n input + 1) parameters, while a feedforward NN model with one hidden layer and one output requires many more parameters: (n input × n neuron + n neuron )+n neuron +1, where n neuron is the number of neurons in the hidden layer. A case of an NN model with varying n neuron will also be investigated. The number of parameters in the model is often used as a proxy for the model complexity.

In the following, the sensitivity of the yield model to the model selection is studied using the LOO and LTO approaches.

Results

Yield model selection -Focus on Bas-Rhin and Landes

In this section, we describe how we selected the most appropriate statistical model for grain maize using the LTO approach based on 22 years of yield data. This test is done over Bas-Rhin and Landes (i.e. two major grain-maize-producing departments in France). We focus here on the LTO results for different models with various selections: number of inputs, model types, number of neurons in the hidden layer, and the number of potential predictors.

We fixed the number of potential predictors n pre = 12 and gradually increased the number of inputs from one to six on the horizontal axis of Fig. 5.2. In both the Bas-Rhin and the Landes examples, underfitting occurs when models are too simple, for example, with only one input. With a higher number of inputs, the LTO procedure shows a similar behaviour as the previous examples from Sect. 4.3.2: the training and validation errors decrease gradually, while the testing errors show an opposite trend. This behaviour suggests that a simple model (e.g. LIN3 for both Bas-Rhin and Landes) is more appropriate.

More complex models were tested in Fig. 5.3: (a) NN models (with 12 potential predictors and 7 neurons in the hidden layer) by increasing the number of inputs and (b) NN3 models (with four potential predictors) by increasing the number of neurons in the hidden layer. The impact of overfitting (Sect. 3.2.2) is noticeable when the model is too complex. For instance, in both cases (Fig. 5.3), the training errors get smaller -close to 0 -for more inputs or more neurons in the hidden layer, as expected. However, the testing and validation errors show large fluctuations when increasing the model complexity. These fluctuations imply that the model is overfitted, and thus, random error or noise appear. Similar results (not shown) were obtained for NN3 models with n potential predictors, where n = 3, 7, 12. Thus, the NN models are unreliable due to the limited number of samples to train a non-linear model. We also tested other examples with LIN3 and NN3 models (Fig. 5.4) to illustrate the cases where model types and the number of potential predictors affect the model quality. Figure 5.4 describes the RMSE values of the predicted grain maize yield anomalies for three datasets (training, validation, and testing) of the LTO procedure. The results of LIN3 models are presented in Fig. 5.4a, and NN3 models (with seven neurons in the hidden layer) are in Fig. 5.4b, with a different number of potential predictors ranging from 3 to 12 on the horizontal axis. The same behaviours are observed: the validation and training errors decrease, while the testing errors increase with the number of potential predictors. Also, the NN3 models show much higher testing and validation RMSE values compared to the LIN3 models. Again, we can conclude -in this grain maize application -that a simpler model will be more beneficial than the complex one.

Reliability model assessment

In this section, a statistical yield model is applied first over 96 French departments to assess the true model quality. Then, we will focus on 10 major departments to assess how the selected models perform for yield anomaly predictions. Figure 5.5 shows the true testing RMSE maps of predicted grain maize yield anomalies in France. Here, the testing errors induced from the LTO procedure are used on the models chosen by the LOO and LTO approaches. In other words, both methods (LOO and LTO) can be considered to identify optimal crop models, but only the LTO method is used (as a reliable tool) to estimate the model generalization ability. For example, when considering only LIN3 models, LOO chooses models with 12 potential predictors, while LTO chooses 3. From these choices, the true model generalization scores (i.e. testing errors) are estimated using the LTO approach, shown in the RMSE maps of Figs. 5.5a1 and b1. Another example focuses on LIN5 models (presented in Figs. 5.5a2 and b2). The true errors obtained from the LOO choice are clearly higher than those from the LTO choice for LIN3 models. For instance, the testing RMSE values range from 10 % to 18 % in many departments in Fig. 5.5a1, while these values are often lower than 10 % in Fig. 5.5b1. This difference shows that the LOO approach underestimates these true errors, as seen in Fig. 5.5a1. Thus, the model choice of the LOO approach is misleading. For more complex models like LIN5 models -that are preferred by the LOO choice -in the second row of Fig. 5.5, higher errors are observed, especially for LOO model errors of many northern departments with up to 22 % of RMSE (Fig. 5.5a2). This grain maize application confirms the benefit of LTO in selecting and assessing the true quality of statistical yield models, while LOO is misleading by underestimating the true errors of its selected models. A simple LIN3 model with three potential predictors is appropriate for this application, considering the limited amount of data.

We now analyse how good the LTO testing estimations are compared to the observations over 10 major grain-maize-producing departments (as shown in Fig. 5.1d). Figure 5.6 presents the boxplots of residuals for these departments, which are the differences between the observed and estimated yield anomalies (Residual = a obs -a est in %). The medians of the residuals lie near zero. This means that the selected models can predict the yield anomalies with acceptable coverage and precision. Although there are some extreme values (Lot-et-Garonne) and some outliers, the interquartile, which ranges from about -8 % to 8 %, shows slight differences between the observations and estimations over study departments.

Seasonal yield forecasting

The LTO approach is helpful for selecting an appropriate model with better forecasting. Here, the model chosen by the LTO procedure is tested for seasonal forecasting, from the sowing time (April) to the forecasting months (i.e. from June to September): all-weather variables (including P and T) from April to June can be selected for the June forecasting. Table 5.1 represents the correlations between the observed and estimated yield anomalies of the forecasts from June to September. The quality of the seasonal forecasting models gradually increases when approaching the harvest because more information is provided. With the weather information at the beginning of the season (April, May, and June), the June forecasting model obtains an average correlation of 0.35 between the observations and estimations. This score is significantly improved when adding information from July (correlation of 0.51). This improvement means that the weather in July strongly influences grain maize yields. The improvement from July to August is much less than from June to July, with an average increase of 0.01 and 0.16, respectively. No information is added in the September forecasting model since it coincides with the harvest time. In other words, the final model should consider only variables from April to August. As in our case, the statistical model selects {T Jul , P M ay , P Apr } as the final inputs for grain maize in the eastern region (Bas-Rhin, Haut-Rhin), {T Jul , P Jul , T Apr } for the southern region (e.g. Landes, Pyrénées-Atlantiques, Gers), and {P Jul , P Apr , P Jun or T Jun } for the central part (Vendée, Charente-Maritime, Vienne). It is reasonable to have different inputs for different regions (or even departments) due to their distinct environmental conditions. In general, weather variables in July -the flowering period -are among the most influential variables. During this time, a high temperature affects the photosynthesis process, thus reducing the potential yield; in contrast, positive precipitation anomalies are preferable [START_REF] Ceglar | Impact of meteorological drivers on regional inter-annual crop yield variability in France[END_REF]Mathieu and Aires, 2018a]. Precipitations in April and May also significantly impact grain maize as water deficit during this vegetative stage decreases plant height [Çakir, 2004].

In addition, Fig. 5.7 shows time series plots of the yield anomaly observations and estimations for different forecasting months in Landes (France). In this case, the June forecasting results show a high correlation with the observed yield anomalies (0.63). This score slightly increases when approaching the harvest. It also indicates that the weather can explain more than 40 % (0.67 2 = 44.89 %) of variations in grain maize yield anomalies in this region, which is in line with other crop studies [START_REF] Ray | Climate variation explains a third of global crop yield variability[END_REF][START_REF] Ceglar | Linking crop yield anomalies to large-scale atmospheric circulation in Europe[END_REF]. However, the forecasting models cannot predict all the extremes (e.g. negative yield anomaly in 1990) that are probably influenced by the climate extremes [Hawkins et al., 2013a;[START_REF] Ceglar | Impact of meteorological drivers on regional inter-annual crop yield variability in France[END_REF]. The statistical models could be improved by adding the indices that focus on extreme weather events.

Conclusions of the chapter

This chapter helps to confirm again that LTO is a general approach. We can easily apply it to other crops, for instance, French grain maize. In detail, the proposed cross-validation technique was used here to choose between simple linear models and more complex neural network models. Our findings also show that a simple LOO is misleading as it overestimates the testing scores. LTO indicated that a simple linear model is preferable because it has a lower testing error. This approach can also be helpful in seasonal forecasting applications (during the growing and the beginning of harvest seasons). In this application, the weather can explain more than 40 % of the variability in yield anomaly, which is similar to that reported in the literature [START_REF] Ray | Climate variation explains a third of global crop yield variability[END_REF][START_REF] Ceglar | Linking crop yield anomalies to large-scale atmospheric circulation in Europe[END_REF]. This score can vary depending on study regions because some regions are more sensitive to the climate than others. Generally, grain maize yield anomalies are mainly influenced by weather variables during the flowering period (July) and the early season (April).

Forthcoming work could develop this LTO into an automatic framework applied to different crops and/or multi-crops. For instance, a multi-crops model on French crops, e.g. spring oat, grain maize, sunflower, sorghum, and soybean could be developed, as they have quite similar growing periods (often ranging from March to October) and data availability for more than 20 years (from 1989 to 2010, at least). This approach will thus result in more samples. In addition, the link between crops can be exploited, for instance, as a regularization technique to reduce the complexity of models used. Additionally, mixed-effects models can be considered instead of a straightforward statistical model. This approach -which intends to use samples in several regions (e.g. gathering samples into groups) to compensate for the lack of historical data -could help us obtain more complex crop models [START_REF] Mathieu | Statistical weather-impact models: An application of neural networks and mixed effects for corn production over the United States[END_REF]. Regarding applications, the crop models we derived here could be used in climate simulations (from an ensemble of climate models for the next 50 years) to investigate the crop yield sensitivity to climate change. In addition to French grain maize, other (annual) crops will be investigated, e.g. over France, Europe, or globally. This task, however, could be very challenging for coffee (both Arabica and Robusta) due to its complex phenology and limited yield database. In the next part of this thesis, we will therefore perform a less ambitious task for coffee suitability. The goal is to estimate future climatically suitable areas for coffee production under several climate change scenarios but not to estimate the future yields.

In the two previous chapters, we considered several statistical impact models to analyse observed yield-climate interactions. The rest of this thesis will focus on future, in particular climate change impacts on coffee production. To accomplish this task, climate simulations are our main ingredient. Nevertheless, climate models are not perfect due to structural or parametric uncertainties, as discussed in Sect. 2.2.2.4. Therefore, we cannot use these climate simulations directly for climate change impact applications. Calibration techniques are introduced to correct these uncertainties and ultimately make the climate simulations fit for application. Many previous studies have applied bias correction or adjustment strategies to remove systematic model errors and calibrate model outputs against observations [START_REF] Wood | Hydrologic implications of dynamical and statistical approaches to downscaling climate model outputs[END_REF][START_REF] Ines | Bias correction of daily gcm rainfall for crop simulation studies[END_REF][START_REF] Lenderink | Estimates of future discharges of the river rhine using two scenario methodologies: direct versus delta approach[END_REF][START_REF] Li | Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change ar4 models using equidistant quantile matching[END_REF][START_REF] Hawkins | Calibration and bias correction of climate projections for crop modelling: An idealised case study over Europe[END_REF][START_REF] Cannon | Bias correction of gcm precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes[END_REF][START_REF] Kim | Bias correction of daily precipitation over South Korea from the long-term reanalysis using a composite Gamma-Pareto distribution approach[END_REF]. This chapter aims to provide a comprehensive reference of available calibration methods for the readers, beginners in the field, as well as the experimented users. Methods will be presented in a pedagogic way -based on the understanding of principles and requirements that lead to mathematical properties -to convey better their ideas, hypotheses, advantages, and inconveniences. Practical diagnostic tools that measure inherent statistical properties are necessary to assess various available methods.

In the following, Sect. 6.1 introduces some useful terminologies, the study area, the datasets, and the notations used. Section 6.2 proposes a synthetic review of the three main calibration "families": direct, delta, and combined. The general principle of each method is described, together with the mathematical formulas, resulting in properties and literature references. In addition, a general table synthesising all this information is provided. Section 6.3 analyses, compares and evaluates six commonly used quantile-based methods. The temporal configuration of the calibration is discussed in Sect. 6.4. Finally,Sect. 6.5

concludes this chapter and discusses further prospects.

This thesis chapter is the subject of a manuscript submitted to Climatic Change journal entitled "Revisiting the bias correction of climate models for impact studies" [Dinh and Aires, 2022b].

Materials

Some useful terminologies

The term "bias" refers to the general concept of "errors" in climate studies. The differences between simulated and observed fields can be characterised among others by the difference in mean, variance, or recurrence interval value [START_REF] Maraun | Towards process-informed bias correction of climate change simulations[END_REF][START_REF] Ivanov | Climate model biases and modification of the climate change signal by intensity-dependent bias correction[END_REF]. The bias can also be defined as the systematic time-independent error component of a model [START_REF] Cannon | Chapter 5 -bias correction of climate model output for impact models[END_REF].

On the other hand, the term "bias correction" is used differently in other fields of research: numerical weather prediction (perfect prognosis and model output statistics), climate model predictions and simulations. In climate studies, several so-called "bias correction" methods were initially developed to correct the biases in the mean and/or variance of climate model outputs. With the same idea of correcting the climate simulations, more complex methods were proposed to modify the entire distribution of the variables to match the observations (i.e. quantilebased techniques), but the term "bias correction" was kept. There is no method available to correct all the biases, however, the quantile-based methods correct more than just the biases in mean and variance but adjust the whole distribution of a given climate variable. An alternative is to use the term "bias adjustment" [START_REF] Addor | Propagation of biases in climate models from the synoptic to the regional scale: Implications for bias adjustment[END_REF][START_REF] Casanueva | Testing bias adjustment methods for regional climate change applications under observational uncertainty and resolution mismatch[END_REF][START_REF] Maraun | Regional climate model biases, their dependence on synoptic circulation biases and the potential for bias adjustment: A process-oriented evaluation of the austrian regional climate projections[END_REF]. In the following, we will refer to "calibration" methods as the general term to transform the climate model simulations towards more realistic values.

Study area

This chapter will focus on the [8 • W -31.5 • E; 33 • N -70

• N] area (shown in Fig. 6.1) that covers Europe and a small part of Africa. It is a very contrasted domain with cold and very hot regions and different elevations. The climate varies a lot from the north to the south, and from the west to the east. For example, the mean temperature difference between the western and eastern parts can be about 10 • C, and this value for monthly precipitation is about 400 mm in January 2014 (shown in Figs. 6. 1b andc). In the Alps, where the elevation is higher than 3 500 m, the January temperature is much cooler than the surrounding area.

Climate data

To account for the seasonality, it is common to have calibration models for each month independently. Therefore, in the following of this chapter, we will consider monthly climate data.

Observations

The re-analysis data, ERA5-Land (Sect. 2.2.1), will be considered in the following as observations. These data are upscaled from their original 0.1 • × 0.1 • resolution to 1 • × 1 • resolution to be consistent with the climate model data. We studied here the gridded monthly-mean 2 m temperature (T) and total precipitation (P) during the 1981-2014 period, which will be referred to in the following as the historical (i.e. present-day) period.

Climate model data

We considered the monthly mean of temperature and precipitation data from CNRM-CM6-1-HR [START_REF] Voldoire | Cnrm-cerfacs cnrm-cm6-1-hr model output prepared for cmip6 highresmip[END_REF][START_REF] Voldoire | Evaluation of cmip6 deck experiments with cnrm-cm6-1[END_REF] as the climate model data (introduced in Sect. 2.2.2.2). They are available from 1850 to 2100 at a resolution of 1 • × 1 • . As mentioned above, the data from 1981-2014 are chosen for the historical period. The model gives the data from 2015 to 2100 as the future projections. These future data are available for several shared socioeconomic pathways (SSPs) [START_REF] O'neill | The Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6[END_REF][START_REF] Riahi | The shared socioeconomic pathways and their energy, land use, and greenhouse gas emissions implications: An overview[END_REF]. In the following, we will perform the tests on the high fossil-fuel development SSP585 scenario (see more in Sect.2.2.2.1). The CNRM model is considered here, but other models could equally be used. This will be done with an ensemble of models in the future, for instance, in chapter 7.

Notations

In the following, the subscript O stands for observation, M for model, C for calibration; h for historical, f for future; and X here refers to T or P variables.

To perform a calibration, three datasets are considered:

• a dataset of historical observations X O,h (i.e. the 1981-2014 period);

• a dataset of historical simulations from the considered climate model X M,h (i.e. the same period as the observations);

• a dataset of model simulations of the future climate (i.e. climate model projections) X M,f (e.g. the 2066-2099 period).

The calibrated values will be denoted as X C,h for values in the historical period and X C,f in the future period.

Distributions

Let us first introduce some basic ideas of several distribution laws that will be used in the following, for different calibration methods.

Normal distribution

The Normal (or Gaussian) distribution can be described by its probability density function (PDF):

f N (x; µ, σ) = 1 σ √ 2π e -1 2 ( x-µ σ ) 2 ,
where µ is the mean or expectation of the distribution, σ is its standard deviation.

The corresponding cumulative distribution function (CDF) is:

F N (x; µ, σ) = 1 2 1 + erf x -µ σ √ 2 ,
where erf is the related error function, e.g. erf (y) = 2 √ pi y 0 e -t 2 dt.

Gamma distribution

The PDF of the Gamma distribution in the shape-scale parametrization is:

f Γ (x; k, θ) = x k-1 e -x θ Γ(k)θ k for x > 0; k, θ > 0,
where Γ(k) is the Gamma function evaluated at k (for all positive integers, Γ(k) = (k -1)!), k is a shape parameter and θ is a scale parameter. The corresponding CDF can be computed as:

F Γ (x; k, θ) = x 0 f Γ (u; k, θ)du = γ k, x θ Γ(k)
where γ k, x θ is the lower incomplete Gamma function.

Generalized Pareto distribution

The generalize Pareto distribution (GPD) is specified by three parameters: location µ, scale σ, and shape ζ [ [START_REF] Dargahi-Noubary | On tail estimation: An improved method[END_REF][START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF]. The PDF of this distribution type is defined by:

f GP D (x; µ, σ, ζ) = 1 σ (1 + ζz) -1 ζ+1
where z = x-µ σ . The corresponding CDF is:

F GP D (x; µ, σ, ζ) = 1 -(1 + ζz) -1 ζ

Calibration of climate simulations 6.2.1 Previous classifications of the calibration methods

In the literature, calibration methods have been classified into distinct categories. [START_REF] Gudmundsson | Technical note: Downscaling rcm precipitation to the station scale using statistical transformations -a comparison of methods[END_REF] defined parametric and non-parametric transformations as the two common calibration methods. [START_REF] Ho | Calibration strategies a source of additional uncertainty in climate change projections[END_REF] determined two calibration pathways: bias correction (i.e. mapping of future climate model simulations to find corresponding calibrated data) and change factor (i.e. mapping of the observations to find future calibrated data). The same classification was identified in later studies, but with different names, e.g. bias correction and delta change in [START_REF] Räisänen | Projections of daily mean temperature variability in the future: cross-validation tests with ensembles regional climate simulations[END_REF], or direct and delta change in [START_REF] Maraun | Bias correcting climate change simulations -a critical review[END_REF]. Another study of [START_REF] Watanabe | Intercomparison of bias-correction methods for monthly temperature and precipitation simulated by multiple climate models[END_REF] used two distinct categories to classify available calibration methods: by whether future statistics are included in the equation for calibration (constant and variable types), and by assumptions on the statistical distributions (parametric and non-parametric types).

Here we divide available methods into three main approaches: direct1 , delta 2 , and combined (as shown in Fig. 6.2). While the direct approach (Sect. 6.2.2.1) applies the calibration directly to the climate model data, the delta approach (Sect. 6.2.2.2) considers only the model change signals (i.e. the change or difference between historic model simulations and future projections). In addition, the combined approach (Sect. 6.2.2.3) associates the two former ones: it calibrates the climate model data but also accounts for the model change signals (e.g. changes in the mean, variance, long-term trend, or changes in quantiles). For example, Fig. 6.2c presents an illustration of a particular combined method (named as equidistant quantile mapping (EqQM), which will be explained in detail in Sect. 6.2.2.3): this method calibrates the raw future model simulations by using the projected model changes in quantiles.

Review of three calibration approaches

Direct approach

The direct approach first estimates the statistical/distributional properties of the climate data (both observations and climate model data) during the historical record (step (1) in Fig. 6.2a). Then it tries to correct the model data to get closer to the observations (step (2) in Fig. 6.2a) to obtain the same statistical properties. The choice of statistical properties depends on the requirements of the considered FIGURE 6.3: Direct calibration methods applied for August temperature, at the pixel [48 • N -7 • E], for the historical and future (2066-2099) periods. LS method: panel (a1) plots the scatter between the observed (X O,h ) and historical model data (X M,h ), and the linear fit line as for X C,h showed in Eq. ( 6.1); panel (a3) shows the PDFs of the observed, modelled and calibrated data (X C,h ); panel (a4) shows the observed data (for comparison purpose), the future model and corresponding calibrated data (X C,f ). EQM method: panel (b1) presents the polynomial fit between the historical model quantiles and the difference between the observed and historical model quantiles; panel (b2) plots the CDFs of the observed, historical model, and calibrated data. From the polynomial fit in panel (b1), we obtain a transfer function to correct the CDF of model data to match the CDF of the observations; panels (b3) and (b4) are the same as panel (a3) and (a4). FQM method: similar to the EQM, but panel (c1) presents a quantile-quantile (Q-Q) plot between the distribution of the observed and historical model data against the expected normal distribution, panel (c2) shows the fitted CDFs instead of the empirical CDFs as in panel (b2).

where µ represents for the historical monthly mean, and i is the time step. The calibrated precipitations are:

P C,h (i) = P M,h (i) × µ (P O,h ) µ (P M,h ) P C,f (i) = P M,f (i) × µ (P O,h ) µ (P M,h ) (6.2)
Since these transformations are additive and/or multiplicative, the term "linear scaling" (LS) is used here. LS corrects the temperature mean, and both precipitation mean and variance (but keeps their ratio constant). This method uses the model data as the baseline; thus, climate variability is more consistent with the model simulations. However, for that reason, the calibrated data is sensitive to the quality of climate model data used as inputs.

Figure 6.3a1 presents the scatter plot between the observed and simulated historical temperatures. Then the PDFs of the observed, modelled, and calibrated temperatures are plotted in Figs. 6.3a3 and a4 for the historical and future periods.

In addition to the LS method, the local intensity scaling method was introduced by [START_REF] Schmidli | Downscaling from gcm precipitation: a benchmark for dynamical and statistical downscaling methods[END_REF] and successfully applied by [START_REF] Moron | Weather types and rainfall over senegal. part ii: Downscaling of gcm simulations[END_REF] or [START_REF] Themeßl | Empirical-statistical downscaling and error correction of daily precipitation from regional climate models[END_REF] to correct the wet-day frequencies and wet-day intensities.

(b) Non-linear scaling

While the LS method focuses on the bias in the mean, other attempts to correct both the mean and variance were considered in the past, for example, using a non-linear transformation model. The non-linearity is not related to the fact that we look at something else than the bias. The power transformation (PT) and modified PT methods [START_REF] Leander | Resampling of regional climate model output for the simulation of extreme river flows[END_REF][START_REF] Leander | Estimated changes in flood quantiles of the river meuse from resampling of regional climate model output[END_REF][START_REF] Van Pelt | Discharge simulations performed with a hydrological model using bias corrected regional climate model input[END_REF][START_REF] Terink | Evaluation of a bias correction method applied to downscaled precipitation and temperature reanalysis data for the rhine basin[END_REF][START_REF] Teutschbein | Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods[END_REF][START_REF] Smitha | An improved bias correction method of daily rainfall data using a sliding window technique for climate change impact assessment[END_REF] consider an exponential form of precipitation by transforming precipitation value (P ) into a corrected amount P c = a × P b . The two parameters a and b in the transformation model make it possible to fit both the mean and the variance.

With a similar idea, Chen et al. [2011a,b] proposed a method aiming at conserving the mean and standard deviation of temperature variables, named as variance scaling. A detailed description of these methods can be found in [START_REF] Teutschbein | Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods[END_REF].

(c) Quantile mapping (QM)

The QM3 is the most widely used technique [START_REF] Panofsky | Some applications of statistics to meteorology[END_REF][START_REF] Wood | Hydrologic implications of dynamical and statistical approaches to downscaling climate model outputs[END_REF][START_REF] Ines | Bias correction of daily gcm rainfall for crop simulation studies[END_REF][START_REF] Déqué | Frequency of precipitation and temperature extremes over France in an anthropogenic scenario: Model results and statistical correction according to observed values[END_REF][START_REF] Boé | Statistical and dynamical downscaling of the seine basin climate for hydro-meteorological studies[END_REF][START_REF] Piani | Statistical bias correction for daily precipitation in regional climate models over Europe[END_REF][START_REF] Watanabe | Intercomparison of bias-correction methods for monthly temperature and precipitation simulated by multiple climate models[END_REF][START_REF] Cannon | Bias correction of gcm precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes[END_REF]. There are numerous QM-based or quantile-based adaptations (as presented in the following sections). The main idea of QM is to correct the distribution function of climate model simulations to match the distribution function of the observations. This is more ambitious than trying to adjust only the mean and variance; we expect here to adjust the whole distribution of the considered variable. A more complex transfer function needs to be used to calibrate data for that purpose. The transfer function will non-linearly adjust all events and thus, be able to correct the frequencies and intensities.

The general QM can be expressed as:

X C,h (i) = F -1 O,h [F M,h (X M,h (i))] X C,f (i) = F -1 O,h [F M,h (X M,f (i))] , (6.3)
where F is the CDF and F -1 is the inverse of F . In Eq. ( 6.3), if F is done without any assumption on the variable distribution (e.g. temperature or precipitation), the method is known as the empirical quantile mapping (EQM) method [START_REF] Themeßl | Empirical-statistical downscaling and error correction of daily precipitation from regional climate models[END_REF]. In this case, F and F -1 are non-parametric and can be quite complex. Figure 6.3 illustrates the EQM method in the second row. Here the transfer function is a polynomial regression that fits the historical model quantiles and the difference between the observed and historical model quantiles (Fig. 6.3b1). This transfer function corrects the model CDF to match the observed CDF during the historical period (Fig. 6.3b2). Then, Fig. 6.3 plots the PDFs of the observed, modelled, and calibrated data in panels (b3) and (b4). By definition, this method is set up during the historical period, and this will conditionate the range of values that can be considered by the calibration model. Thus, the problem arises when a simulated value is out-of-range of the historical values leading to saturation effects. For instance, the fit on this overrange data point can cause "new extremes". To deal with these extremes, [START_REF] Wood | Hydrologic implications of dynamical and statistical approaches to downscaling climate model outputs[END_REF] and [START_REF] Themeßl | Empirical-statistical downscaling and error correction of regional climate models and its impact on the climate change signal[END_REF] proposed some form of extrapolation using a parametric distribution based on a hypothesis for the variable distribution. [START_REF] Boé | Statistical and dynamical downscaling of the seine basin climate for hydro-meteorological studies[END_REF] employed the constant correction approach to handle this new extreme issue. In another study, [START_REF] Michelangeli | Probabilistic downscaling approaches: Application to wind cumulative distribution functions[END_REF] suggested a specific extreme value method to model extreme wind intensities.

Instead of using the non-parametric CDF (i.e. EQM), another suggestion is to fit the CDF to the common distribution laws, which we call here a fitted quantile mapping (FQM) method. This method performs best when the distribution of the data fits well the distribution type. For example, the Normal distribution N (µ, σ) is often used to fit the temperature distribution [START_REF] Thom | The rational relationship between heating degree days and temperature[END_REF][START_REF] Hay | Use of statistically and dynamically downscaled atmospheric model output for hydrologic simulations in three mountainous basins in the western united states[END_REF][START_REF] Haerter | Climate model bias correction and the role of timescales[END_REF][START_REF] Ho | Calibration strategies a source of additional uncertainty in climate change projections[END_REF][START_REF] Teutschbein | Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods[END_REF].

On the other hand, the two-parameter gamma distribution [START_REF] Thom | A Note on the Gamma Distribution[END_REF] Γ(k, θ) has been shown to be appropriate for precipitation distribution [START_REF] Watterson | Simulated changes due to global warming in daily precipitation means and extremes and their interpretation using the gamma distribution[END_REF][START_REF] Hay | Use of statistically and dynamically downscaled atmospheric model output for hydrologic simulations in three mountainous basins in the western united states[END_REF][START_REF] Ines | Bias correction of daily gcm rainfall for crop simulation studies[END_REF][START_REF] Block | A streamflow forecasting framework using multiple climate and hydrological models1[END_REF][START_REF] Piani | Statistical bias correction for daily precipitation in regional climate models over Europe[END_REF].

For temperature, the calibrated model data are given by:

T C,h (i) = F -1 N F N (T M,h (i); µ (T M,h ) , σ (T M,h ) ); µ (T O,h ) , σ (T O,h ) T C,f (i) = F -1 N F N (T M,f (i); µ (T M,h ) , σ (T M,h ) ); µ (T O,h ) , σ (T O,h ) (6.4)
where F N is the Normal CDF and F -1 N is its inverse. For precipitation, we use:

P C,h (i) = F -1 Γ F Γ (P M,h (i); k (P M,h ) , θ (P M,h ) ); k (P O,h ) , θ (P O,h ) P C,f (i) = F -1 Γ F Γ (P M,f (i); k (P M,h ) , θ (P M,h ) ); k (P O,h ) , θ (P O,h ) (6.5)
where F Γ is the gamma CDF and its inverse is denoted as F -1 Γ . The methodology of the FQM method is also presented in the third row of Fig. 6.3: panel (c1) presents a quantile-quantile (Q-Q) plot between the distribution of the observed and model data during the historical record against the expected normal distribution. After that, the fitted model CDF is transferred to the observed CDF in panel (c2). Finally, the PDFs of all data are plotted for the historical and future periods in panels (c3) and (c4).

Other distribution laws (e.g. Gumbel, Weibull, Lognormal, etc.) can also be considered to better fit the variables of interest. For example, the generalized extreme value (GEV) distribution can be used to fit the temperature and precipitation extremes [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF][START_REF] Kharin | Changes in temperature and precipitation extremes in the ipcc ensemble of global coupled model simulations[END_REF]. [START_REF] Watterson | Calculation of probability density functions for temperature and precipitation change under global warming[END_REF] showed that the four-parameter beta distribution provides a smooth PDF matching the mean and range of the simulated data.

The Gamma-Pareto quantile mapping -based on the combination of a gamma distribution and a generalized Pareto distribution GP D(k, θ, ζ) [START_REF] Dargahi-Noubary | On tail estimation: An improved method[END_REF][START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF] -is widely used in the literature [START_REF] Gutjahr | Comparing precipitation bias correction methods for high-resolution regional climate simulations using cosmo-clm[END_REF][START_REF] Volosciuk | A combined statistical bias correction and stochastic downscaling method for precipitation[END_REF][START_REF] Kim | Bias correction of daily precipitation over South Korea from the long-term reanalysis using a composite Gamma-Pareto distribution approach[END_REF]. The main idea of this method is to better correct the higher percentiles. For precipitation, F Γ is fitted for values between the 5 th and 95 th percentiles, while the GPD is used to fit the upper and lower 5 %:

P C,f (i) = F -1 Γ F Γ (P M,f (i); k (P M,h ) , θ (P M,h ) ); k (P O,h ) , θ (P O,h ) , if 5 th percentile < P M,f (i) < 95 th percentile ; P C,f (i) = F -1 GP D F GP D (P M,f (i); µ (P M,h ) , σ (P M,h ) , ζ (P M,h ) ); µ (P O,h ) , σ (P O,h ) , ζ (P O,h ) , if P M,f (i) ≥ 95 th percentile or P M,f (i) ≤ 5 th percentile. (6.6)
Here, the GPD is specified by three parameters: location µ, scale σ, and shape ζ, as presented in Sec. 6.1.5.

The formula of Eq. ( 6.6) can also be used to calibrate the temperature variable, but a Normal distribution F N will be preferred over F Γ for the 5 th and 95 th percentiles.

We have shown that a different transfer function needs to be used in the direct approach, depending on the statistical properties that we aim to preserve in the calibrated data (e.g. mean, variance, quantiles, extremes).

Delta approach

The delta approach extracts the model change signals, i.e. the differences between historical and future climate simulations. These signals are then transferred into the observations to obtain the calibrated climate projections. Thus, this delta approach is not a calibration or bias correction of climate models [START_REF] Teutschbein | Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods[END_REF][START_REF] Maraun | Bias correcting climate change simulations -a critical review[END_REF]. However, as it has a long history and has been widely used in climate studies [START_REF] Hay | A comparison of delta change and downscaled gcm scenarios for three mountain basins in the united states[END_REF]Graham et al., 2007a,b;[START_REF] Olsson | Applying climate model precipitation scenarios for urban hydrological assessment: A case study in kalmar city, sweden[END_REF][START_REF] Bosshard | Spectral representation of the annual cycle in the climate change signal[END_REF][START_REF] Hawkins | Calibration and bias correction of climate projections for crop modelling: An idealised case study over Europe[END_REF][START_REF] Navarro-Racines | High-resolution and bias-corrected cmip5 projections for climate change impact assessments[END_REF], we will therefore include it in this review study as one of the available climate calibration approaches.

Like the direct approach (Sect. 6.2.2.1), the delta approach can be used to adjust different statistical properties (e.g. mean, variance, or quantiles). All methods presented in Sect. 6.2.2.1 can be adapted to the delta approach in practice. Here, we present the two main delta methods: LS ∆ and QM ∆ (the subindex ∆ will be added to avoid confusion with the direct approach). In addition, an adjustment of QM -named the adjusted quantile mapping (AQM) -will be introduced.

Since, by definition, this approach does not aim to correct the historical model data, only the calibrations for the future period X C,f will be presented in the following of this section.

(a) Linear scaling (LS ∆ )

In the LS ∆ method [Graham et al., 2007a;Chen et al., 2011a;[START_REF] Hawkins | Calibration and bias correction of climate projections for crop modelling: An idealised case study over Europe[END_REF], the absolute (or relative) changes in temperature (or precipitation) characteristics, estimated from climate model data, are transferred to the observed temperature (or precipitation) time series by additive (or multiplicative) factors:

T C,f (i) = T O,h (i) + µ (T M,f ) -µ (T M,h ) (6.7)
for temperature, and:

P C,f (i) = P O,h (i) × µ (P M,f ) µ (P M,h ) (6.8)
for precipitation. Since this method uses observed climate as a baseline, the obtained calibrations will conserve all details of the observations. For instance, the number of rainy days will not change for future climate. In addition, extreme events are modified by the same factors as all other precipitation events.

(b) Quantile mapping (QM ∆ )

The QM ∆ method assumes that the changes from historical to future in the observation distribution are the same as the changes in the model distribution [START_REF] Olsson | Applying climate model precipitation scenarios for urban hydrological assessment: A case study in kalmar city, sweden[END_REF][START_REF] Ho | Calibration strategies a source of additional uncertainty in climate change projections[END_REF][START_REF] Hawkins | Calibration and bias correction of climate projections for crop modelling: An idealised case study over Europe[END_REF]. This is achieved by using QM (Sect. 6.2.2.1c): the transfer function is first estimated from the changes in the model distribution (between future and historical data); then, this transfer function is applied to the observed time series. While the LS ∆ method corrects only temperature mean and applies the same factors for precipitation variable, the QM ∆ method is able to correct both frequencies and intensities.

In general, the calibrated data (both precipitation and temperature) can be estimated as:

X C,f (i) = F -1 M,f [F M,h (X O,h (i))] (6.9)
As for the direct approach, the CDF function F can be empirical (EQM ∆ ) or follow a known distribution law (FQM ∆ ) (see Sect. 6.2.2.1).

(c) Adjusted quantile mapping (AQM)

As an adaptation of the QM ∆ method, [START_REF] Amengual | A statistical adjustment of regional climate model outputs to local scales: Application to platja de palma, spain[END_REF] proposed the AQM method (or quantile-quantile adjustment) to obtain future model simulations. The method first detects the changes of each quantile in the CDFs of the climate model outputs and then applies these changes, after rescheduling based on the historical period, to the observed time series.

The mathematical formula of the calibrated future model data is as follows:

X j C,f = X j O,h + gδ + s(δ j -δ) (6.10)
where:

• the index j here refers to the j th ranked value of the corresponding CDFs (e.g. for observations or historical and future model data),

•

δ j = X j M,f -X j M,h , • δ = N j=1 δ j N , • g = 1 for temperature and g = µ (X O,h ) µ (X M,h ) for precipitation, • s = IQR| O,h
IQR| M,h . Here, IQR| O,h and IQR| M,h are the interquartile ranges of the observed and modelled data during the historical record, respectively.

Combined approach

Various quantile-based methods are based on the QM technique (direct approach in Sect. 6.2.2.1) but also account for the model change signals (e.g. delta approach). This is intended to benefit from the advantages of the two previously described approaches (Sects. 6.2.2.1 and 6.2.2.2). These methods are called here "combined approach".

In this approach, the calibrated model simulations over the historical record are identical to the results obtained by the direct QM method (i.e. X C,h in Eq. ( 6.3)), except for the scaled distribution mapping (SDM) method. Thus, only calibrated model simulations over the future record (i.e. calibrated model projections) will be considered in the following.

(a) New quantile mapping (NewQM)

Suggested by [START_REF] Watanabe | Intercomparison of bias-correction methods for monthly temperature and precipitation simulated by multiple climate models[END_REF], a NewQM method -which is based on FQM -takes into account the changes in mean and standard deviation (or coefficient of variation) for temperature (or precipitation).

Mathematically, the model temperatures over the future period are calibrated following:

T C,f (i) = F -1 N F N (T M,f (i); µ (T M,h ) , σ (T M,h
) ); µ cor , σ cor (6.11) in which

µ cor = µ (T O,h ) + µ (T M,f ) -µ (T M,h ) and σ cor = σ (T O,h ) ×σ (T M,f ) σ (T M,h )
.

The non-precipitation months (P = 0) are first removed, and the model precipitations over the future period are then corrected following:

P C,f (i) = F -1 Γ F Γ (P M,f (i); k (P M,h ) , θ (P M,h ) ); k cor , θ cor (6.12)
where k cor and θ cor are induced from the corrected mean and coefficient of variation (CV), in which:

µ cor = µ (P O,h ) ×µ (P M,f ) µ (P M,h )
and CV cor =

CV (P O,h ) ×CV (P M,f ) CV (P M,h )
.

(b) Detrended quantile mapping (DetQM)

Based on the idea that trends should be preserved during the calibration process [START_REF] Hempel | A trend-preserving bias correction -the isi-mip approach[END_REF], the DetQM was suggested by [START_REF] Cannon | Bias correction of gcm precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes[END_REF]. This method is developed from the EQM method (Sect. 6.2.2.1), but it explicitly accounts for the modelled changes in trend. In detail, the model data are used to obtain a climate trend. Then, QM is used to calibrate model projections that have their long-term trend removed. After that, the removed climate trend is then reintroduced. This DetQM method is designed to preserve the long-term absolute (or relative) trend of the model temperature (or precipitation) data. This method follows three steps:

• Step (1): Trend removal: T Det1 C,f (i) = T M,f (i)-µ (T M,f ) -µ (T M,h ) ; and P Det1 C,f (i) = P M,f (i) × µ (P M,h ) µ (P M,f ) . • Step (2): Quantile mapping: T Det2 C,f (i) = F -1 O,h F M,h (T Det1 C,f (i)) ; and P Det2 C,f (i) = F -1 O,h F M,h (P Det1 C,f (i)) . • Step (3): Trend reimposition: T C,f (i) = T Det2 C,f (i) + µ (T M,f ) -µ (T M,h ) ; and P C,f (i) = P Det2 C,f (i) × µ (P M,f ) µ (P M,h ) .
Here, zero values in the observed and modelled data are first replaced by nonzero uniform random values below a tracing threshold prior to calibration. Then, values -after the calibration -that are lower than the predefined threshold are set back to zero. This measure, therefore, can adjust the wet-day frequencies.

Initially, EQM is considered in step ( 2), but we can also apply FQM instead, depending on the considered application.

(c) Equidistant quantile mapping (EqQM)

The EqQM method [START_REF] Li | Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change ar4 models using equidistant quantile matching[END_REF][START_REF] Pierce | Improved bias correction techniques for hydrological simulations of climate change[END_REF] focuses on the historical bias, i.e. the difference between observed and modelled data, at a given percentile during the historical period will apply to the future period:

X C,f (i) = X M,f (i) + F -1 O,h [F M,f (X M,f (i))] -F -1 M,h [F M,f (X M,f (i))] . (6.13)
Initially, [START_REF] Li | Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change ar4 models using equidistant quantile matching[END_REF] used EqQM to correct both temperature and precipitation data. However, the method was shown to be inefficient for precipitation; for example, it was problematic over some dry regions [START_REF] Li | Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change ar4 models using equidistant quantile matching[END_REF] or resulted in negative precipitation values [START_REF] Wang | Equiratio cumulative distribution function matching as an improvement to the equidistant approach in bias correction of precipitation[END_REF]. Consequently, a multiplicative factor (i.e. ratio) should be considered instead of an additive form, as often for precipitation. This adaptation was first proposed by [START_REF] Li | Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change ar4 models using equidistant quantile matching[END_REF] and then introduced as the equiratio CDF matching method by [START_REF] Wang | Equiratio cumulative distribution function matching as an improvement to the equidistant approach in bias correction of precipitation[END_REF]: .14) For simplicity reasons, we will consider Eqs. (6.13) and (6.14) as the additive and multiplicative forms of the EqQM method.

P C,f (i) = P M,f (i) × F -1 O,h [F M,f (P M,f (i))] F -1 M,h [F M,f (P M,f (i))] . ( 6

(d) Quantile delta mapping (QDM)

The QDM method [START_REF] Cannon | Bias correction of gcm precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes[END_REF] aims to preserve the change signals of the model in terms of quantiles. The method considers the relative changes of a ratio variable. QDM follows three main steps:

• Step (1): The future model simulations are detrended by quantile and adjusted to observations by QM:

P QDM 1 C,f (i) = F -1 O,h [F M,f (P M,f (i))] (6.15)
• Step (2): The model changes (relative for precipitation) in terms of quantiles are computed following:

∆ m,P (i) = F -1 M,f [F M,f (P M,f (i))] F -1 M,h [F M,f (P M,f (i))] = P M,f (i) F -1 M,h [F M,f (P M,f (i))] (6.16)
•

Step (3): The relative changes ∆ m,P (i) are then superimposed on the corrected data

P QDM 1 C,f (i)
:

P C,f (i) = P QDM 1 C,f (i) × ∆ m,P (i) (6.17)
In short, Eq. ( 6.17) can be rewritten as

P C,f (i) = P M,f (i) × F -1 O,h [FM,f(PM,f(i))] F -1 M,h [FM,f(PM,f(i))]
,

which is similar to the equiratio CDF matching method [START_REF] Wang | Equiratio cumulative distribution function matching as an improvement to the equidistant approach in bias correction of precipitation[END_REF] showed in Eq. (6.14). Also, the additive form of the QDM method had been shown to be equivalent to the EqQM as in Eq. ( 6.13) [START_REF] Li | Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change ar4 models using equidistant quantile matching[END_REF][START_REF] Cannon | Bias correction of gcm precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes[END_REF].

Although EqQM and QDM start with different motivations (i.e. the historical biases and the future change signals respectively), they are found to be similar, as concluded in [START_REF] Cannon | Bias correction of gcm precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes[END_REF]. In the following sections, we will consider these two methods (EqQM and QDM) as EqQM.

In addition, QDM uses the same correction for the wet-day frequencies as the DetQM method (shown previously in Sect. 6.2.2.3b).

(e) Scaled distribution mapping (SDM)

The SDM method was first proposed by [START_REF] Switanek | Scaled distribution mapping: a bias correction method that preserves raw climate model projected changes[END_REF] and recently used to correct a regional climate model by [START_REF] Maraun | Regional climate model biases, their dependence on synoptic circulation biases and the potential for bias adjustment: A process-oriented evaluation of the austrian regional climate projections[END_REF]. This method is similar to QDM but also considers the frequency of rain days (which is essential when dealing with the daily precipitation data) and the likelihood of individual events.

SDM can be implemented by following seven steps for both temperature and precipitation:

• Step (1): preprocess data: detrend the raw modelled and observed temperature time series, and set a precipitation threshold to separate days with and without rain;

• Step (2): fit a PDF on the preprocessed data;

• Step (3): calculate the scaling between the fitted raw future model distribution and the fitted raw historical distribution at each probability corresponding to the events of the raw future model time series;

• Step (4): calculate the recurrence intervals (RI);

• Step (5): find the scaled or adjusted RI (RI scaled ) and then the corresponding scaled CDF values (F scaled ) for the future simulation;

• Step ( 6): calculate the initial array of corrected values (BC initial );

• Step (7): reinsert BC initial back into the right time series, and then correct the precipitation and trend of temperature.

A detailed algorithm of SDM can be found in [START_REF] Switanek | Scaled distribution mapping: a bias correction method that preserves raw climate model projected changes[END_REF].

The BC initial can be mathematically written as:

BC initial (i) = F -1 O,h [F scaled,T (T M,f (i))] + ∆ m,T (i) for temperature, (6.18) in which ∆ m,T (i) = F -1 M,f [F M,f (T M,f (i))] -F -1 M,h [F M,f (T M,f (i))] × σ (T O,h ) σ (T M,h
) . The multiplicative form of SDM for precipitation is given by:

BC initial (i) = F -1 O,h [F scaled,P (P M,f (i))] × ∆ m,P (i), (6.19)
in which ∆ m,P (i) is given in Eq. (6.16). As we can see by comparing Eqs. (6.17) and (6.19), the two methods (QDM and SDM) are very similar, except that SDM considers F scaled instead of the CDF of the model projections F M,f . The F scaled reflects the scaling of the model change in likelihood event corresponding to the observed likelihoods. In addition, SDM includes the preprocessing step that also helps dealing with the trend in the temperature data or the non-precipitation day.

As a result, SDM may work better when working with daily precipitation data.

General comments

Table 6.1 synthesises the methods discussed in Sect. 6.2.2, providing the main ideas, their advantages and disadvantages, and references. Each method is based on different assumptions and aims to correct different statistical properties of climate model data (except the delta approach), e.g. simply correcting the mean (LS), the mean and variance (non-linear scaling), or adjusting all moments of the climate model distribution (QM and its variants). Therefore, the choice of using a particular calibration method depends on the applications [START_REF] Ho | Calibration strategies a source of additional uncertainty in climate change projections[END_REF][START_REF] Teutschbein | Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods[END_REF][START_REF] Switanek | Scaled distribution mapping: a bias correction method that preserves raw climate model projected changes[END_REF] and on what the impact modellers aim to emphasise on the observations (e.g. average mean, high-resolution pattern) and on the model data (e.g. trend, distributions of changes).

METHODS

MAIN IDEAS & COMMENTS REFERENCES

DIRECT APPROACH

Applies the calibration to the model simulations; follows the stationary assumption (i.e. the model discrepancies are time-invariant).

Linear Scaling (LS)

LS

LS adjusts the mean values by applying to the raw model data an additive (or multiplicative) term, corresponding to the difference between mean observed and historical modelled data. [START_REF] Lenderink | Estimates of future discharges of the river rhine using two scenario methodologies: direct versus delta approach[END_REF]Graham et al., 2007a;[START_REF] Berg | Bias correction of high resolution regional climate model data[END_REF] LS corrects the temperature mean and both precipitation mean and variance (but keeps their ratio constant). This method offers calibrated data with a variability more consistent with the raw model data. But it adjusts all events by the same correction factor and thus cannot correct the frequencies.

Local intensity scaling

This method has the same idea as LS, but also considers the wet-day frequencies and intensities. [START_REF] Schmidli | Downscaling from gcm precipitation: a benchmark for dynamical and statistical downscaling methods[END_REF][START_REF] Moron | Weather types and rainfall over senegal. part ii: Downscaling of gcm simulations[END_REF][START_REF] Themeßl | Empirical-statistical downscaling and error correction of daily precipitation from regional climate models[END_REF] Non-Linear Scaling Power Transformation (PT); Modified PT PT and modified PT use an exponential form to adjust the variance statistics of precipitation time series. [START_REF] Leander | Resampling of regional climate model output for the simulation of extreme river flows[END_REF][START_REF] Leander | Estimated changes in flood quantiles of the river meuse from resampling of regional climate model output[END_REF][START_REF] Van Pelt | Discharge simulations performed with a hydrological model using bias corrected regional climate model input[END_REF][START_REF] Terink | Evaluation of a bias correction method applied to downscaled precipitation and temperature reanalysis data for the rhine basin[END_REF]Smitha et al., 2018 Variance Scaling This method first corrects the mean, as for LS, then accounts for the variance correction of temperature time series. Chen et al. 2011a,b;[START_REF] Teutschbein | Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods[END_REF] Non-Linear Scaling method can correct both the mean and variance of the raw model data but not the wet-day frequencies and intensities.

Quantile Mapping (QM)

The main idea of QM is to correct the climate modelled distribution function to match the observed distribution function by using a transfer function. Empirical Quantile Mapping (EQM)

The transfer function is done without any assumption of the variable distribution (i.e. nonparametric method). To improve the performance of EQM on its ability to generate new extremes, some studies suggested some form of extrapolations or a constant correction approach. Panofsky and Brier 1986;[START_REF] Wood | Hydrologic implications of dynamical and statistical approaches to downscaling climate model outputs[END_REF][START_REF] Déqué | Frequency of precipitation and temperature extremes over France in an anthropogenic scenario: Model results and statistical correction according to observed values[END_REF][START_REF] Boé | Statistical and dynamical downscaling of the seine basin climate for hydro-meteorological studies[END_REF][START_REF] Michelangeli | Probabilistic downscaling approaches: Application to wind cumulative distribution functions[END_REF]Themeßl et al., 2012 Fitted Quantile Mapping (FQM) FQM try to fit the CDFs to the common known distribution laws (e.g. normal, gamma, GEV, gamma-pareto, or beta). [START_REF] Watterson | Simulated changes due to global warming in daily precipitation means and extremes and their interpretation using the gamma distribution[END_REF][START_REF] Ines | Bias correction of daily gcm rainfall for crop simulation studies[END_REF][START_REF] Block | A streamflow forecasting framework using multiple climate and hydrological models1[END_REF][START_REF] Piani | Statistical bias correction for daily precipitation in regional climate models over Europe[END_REF][START_REF] Ho | Calibration strategies a source of additional uncertainty in climate change projections[END_REF] QM non-linearly adjusts all events and thus, be able to correct the frequencies and intensities. EQM may raise the problem of functional robustness, e.g. causing the new extremes, since its overlarge degree of freedom is equal to the number of parameters. FQM performs best when the data distribution fit well the chosen distribution law. Otherwise, unrealistic values can be produced.

DELTA APPROACH

Assumes that the observed change from historical/present-day to future (e.g. change in mean or distribution) will be similar to the model change. In detail, this approach employs the model response to climate change to modify observations; therefore, it is not a calibration of a climate model.

Linear Scaling (LSΔ )

LSΔ calibrates the future model data by adding (or multiplying) the model change signals (mean bias or ratio) to the observed time series. Graham et al., 2007a;Chen et al., 2011a;Hawkins et al., 2013 Quantile Mapping (QMΔ) QMΔ assumes that the change from historical to future in the observation distribution is the same as the change in the model distribution. [START_REF] Olsson | Applying climate model precipitation scenarios for urban hydrological assessment: A case study in kalmar city, sweden[END_REF][START_REF] Ho | Calibration strategies a source of additional uncertainty in climate change projections[END_REF][START_REF] Hawkins | Calibration and bias correction of climate projections for crop modelling: An idealised case study over Europe[END_REF] Adjusted Quantile Mapping (AQM) AQM first detects the changes of each quantile in the CDFs of model data and then applies these changes, after reschedules based on the historical period, to the observed time series. [START_REF] Amengual | A statistical adjustment of regional climate model outputs to local scales: Application to platja de palma, spain[END_REF] By definition, delta methods have the advantages of conserving all small details of the observations and the change in variabilities of the climate model. However, the method artificially accounts for all observed variabilities that, in fact, should not be resolved in the climate model data.

COMBINED APPROACH

Aims to calibrate directly the climate model data but also accounts for the model change signals; is not constrained by the stationarity assumption (except DetQM).

New Quantile Mapping (NewQM)

NewQM corrects the mean and std (or mean and coefficient of variation), instead of using directly the ones from observations, prior to QM. [START_REF] Watanabe | Intercomparison of bias-correction methods for monthly temperature and precipitation simulated by multiple climate models[END_REF] This method conserves the changes in several statistical parameters from the historical to the future period between uncalibrated and calibrated data.

Detrended Quantile Mapping (DetQM)

DetQM first removes (or rescales) the modelled trend in the long-term mean before QM and then reimposes it afterwards.

Comparison of quantile-based calibration methods

The delta methods -especially LS ∆ and QM ∆ -do not adjust the climate model outputs. However, among three delta methods listed in Sect. 6.2.2.2, it can be interesting to analyse the AQM method, which is based on QM ∆ but also accounts for the relationship between observations and model simulations in the historical period (i.e. the factors g and s in Eq. ( 6.10)). In the following sections, our analyses will focus on six quantile-based methods (i.e. QM and its variants): one direct method (FQM), one delta method (AQM), and four combined methods (NewQM, DetQM, EqQM, SDM). Here, the QM refers to the FQM, where a normal distribution is used to fit temperature values, and a gamma distribution is for precipitation values. Other distribution functions or EQM could have been considered. For the detailed analysis among direct methods (and LS ∆ ), we refer the readers to the previous review study by [START_REF] Teutschbein | Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods[END_REF].

By definition, the historical calibrated data of the combined methods (except for SDM) will be similar to the one from FQM, i.e. X C,h (FQM). Also, to be consistent, X C,h (AQM) will be considered equivalent to the observations X O,h . In brief, the results include X C,h (FQM), X C,h (AQM), and X C,h (SDM) for historical period and six calibration methods for the future period.

Metrics of comparison

The calibration for the future period, in theory, cannot be validated. Previous studies often tried to evaluate the calibration performance on the historical periods, but this evaluation might not be valid for the future [START_REF] Teutschbein | Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods[END_REF]. [START_REF] Switanek | Scaled distribution mapping: a bias correction method that preserves raw climate model projected changes[END_REF] tested the impact of different historical periods (i.e. control periods) and showed that the control period strongly influences the calibration. This is equivalent in machine learning; they try to learn too much on the training data at the detriment to the generalization ability of the model (i.e. overtraining or overfitting). As a result, this study will not evaluate calibration methods based on the historical period. Instead, we focus on the two main objectives (shown in Fig. 6.4) that the calibration is designed for: (1) Do the calibrated data (X C,h ) get closer to the observed data (X O,h )? ( 2 • standard deviation or std (σ) -an indication of variability,

• 10 th and 90 th percentiles (X10, X90) -an indication of extreme values,

• skewness (γ) -a measurement of the asymmetry of the distribution of values.

We also included the trend for temperature analysis and coefficient of variation (CV) for precipitation analysis. Furthermore, two diagnostics are considered: the spatial correlation (COR) and root mean square error (RMSE) between D M od and D Cal for all grid cells. Here we can see that the raw climate model data (X M,h and X M,f ) change significantly over time. For example, the July temperature variability increases about 2 • C over the study area (shown in Fig. 6.5a). In addition, the skewness varies a lot between the two considered periods. In Fig. 6.5b, the skewness of the July temperature series over France is positive for the historical period (panel of X M,h ), i.e. right skewness -most values are on the left side of the mean value. However, this value is mostly negative for the future period (panel of X M,f ). These changes, for instance, in standard deviation and skewness, suggest that a calibration method should better correct higher-order characteristics of the climate distribution rather than just correct the mean [START_REF] Li | Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change ar4 models using equidistant quantile matching[END_REF].

Performance

Impact of calibration methods on climate model outputs

For the considered examples in Fig. 6.5, three quantile-based methods (FQM, AQM, and SDM) are able to correct the standard deviation. For the historical period, FQM calibrates X M,h to be similar to X O,h (i.e. the correlation between X C,h (F QM ) and X O,h is 1). AQM, by definition, considers the observations as the baseline data. SDM reduces the variability but not the same as X O,h . For the future period, all methods are able to reduce the variability. AQM shows very large variability over the northeast regions because of the increase in variability between X M,f and X M,h ; this increase, by definition, is introduced into the observed variability. SDM does improve the model but not so much due to its target that focuses on preserving the raw climate model change (more details will be discussed in the next parts). On the other hand, each method deals with other characteristics of the distribution (e.g. skewness) differently. In Fig. 6.5b, X C,h (F QM ) is exactly the same as X M,h ; and X C,f (F QM ) is quite similar to X M,f . Although FQM adjusts all events, this method underlies an assumption that the climate distribution does not change much over time: it is stationary in the skew of the distribution and only changes in the mean and standard deviation. For the delta method, the skewness of X C,f (AQM ) tends to retain signals/patterns from X O,h , but also accounts for the model change signal. SDM adjusts the climate model data for both historical and future periods. While Figs. 6.6 and 6.7 show particular examples in July, Fig. 6.8 summarises the temporal evolution of six calibration methods over all pixels shown in Fig. 6.6, for both temperature and precipitation. Here, we used a calibration model for each month of the year. For temperature, AQM and four combined methods (NewQM, DetQM, EqQM, SDM) outperform the traditional FQM method in preserving the climate change signals for all considered characteristics and the trend (as seen in Fig. 6.6). For example, FQM shows the RMSE(D M od ,D Cal ) in the change to the mean of about 1.5 to 2 • C while other methods much better minimise this RMSE value (approximately 0 • C). Similar behaviours are observed in other characteristics of the distribution (standard deviation, 10 th and 90 th percentiles) in which FQM shows much higher RMSE values than other methods. SDM performs best with respect to preserving the changes in the trend. This can be explained by the preprocessing step of SDM to adjust the trend of temperature. For precipitation, five quantile-based variants (AQM, NewQM, DetQM, EqQM, and SDM) show improvement in the ability to preserve the raw projected mean change compared to FQM. On the other hand, FQM preserves well the climate change signals for other moments of the distribution. Again, NewQM performs better for the coefficient of variation than standard deviation as its target is to focus on the coefficient of variation correction. There is no significant difference among considered methods for precipitation. This behaviour can be explained by the considered temporal resolution: most methods attempt to handle the non-rain days, while this study focuses on the monthly data.

Impact of calibration methods on climate change signals

To sum up, Fig. 6.8 gives an overview/intercomparison of all methods via several diagnostics. More methods/diagnostics -that are important to the users -can be added easily. This kind of presentation is something simple but very informative. For instance, it is obvious to observe the seasonal patterns of each method or to identify which method is better in preserving the climate change signals (e.g. change in the mean or trend). With this synthetic view, one can automatically choose the suitable method for their applications. We thus suggest that this comparison figure (i.e. Fig. 6.8) should become a standard for assessing calibration methods in climate impact studies.

Spatio-temporal configurations of the calibration

The setting of the calibration model can be done on several space and time configurations. This will impact the properties of the calibration results.

Temporal configuration

Seasonal calibration

For monthly data, a seasonal calibration can be considered. Instead of using the same calibration model for all months, a calibration model is set up here for each month of the year. This is done by compositing, for instance, all the January months to set up a January calibration. This approach is beneficial when the expected corrections on climate simulations follow a seasonal pattern, as the corrections might differ from the summer and winter months. This approach can have several limitations: (1) the sample size is reduced (e.g. only 30 monthly samples for the 1981-2010 period), (2) a large artificial jump can be observed between two consecutive months since each month has its own independent transfer function (but this could be reduced by using a 3-month window), and (3) the precise temporal pattern of the model can be distorted. We thus apply a 3-month moving window [START_REF] Cannon | Bias correction of gcm precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes[END_REF] for our monthly calibration: Fig. 6.9 presents the performance of six calibration methods for all the months from January to December, by using such a temporal moving window instead of the monthly calibration shown in Fig. 6.8. As expected, the transitions appear to be smoother with this approach. In addition, the moving window technique shows a significant improvement in the ability to preserve the climate change signals for both temperature and precipitation, for the six compared methods, as the RMSEs between D Cal and D M od are much lower than the ones in Fig. 6.8. For daily data, more options are available. For example, in [START_REF] Wood | Long-range experimental hydrologic forecasting for the eastern united states[END_REF]; [START_REF] Maurer | The utility of daily large-scale climate data in the assessment of climate change impacts on daily streamflow in california[END_REF], all days in a particular month are composited to set up a calibration model for that particular month. As with monthly data, there is still the problem of discontinuities at the edges of the time window. For instance, the January 31 st is calibrated using the January model and the February 1 st is calibrated using the February model [START_REF] Thrasher | Technical note: Bias correcting climate model simulated daily temperature extremes with quantile mapping[END_REF][START_REF] Pierce | Improved bias correction techniques for hydrological simulations of climate change[END_REF]. A moving window approach, using ± n days instead of a static month, was introduced to avoid this issue [START_REF] Themeßl | Empirical-statistical downscaling and error correction of daily precipitation from regional climate models[END_REF][START_REF] Thrasher | Technical note: Bias correcting climate model simulated daily temperature extremes with quantile mapping[END_REF][START_REF] Gennaretti | Toward daily climate scenarios for canadian arctic coastal zones with more realistic temperatureprecipitation interdependence[END_REF][START_REF] Pierce | Improved bias correction techniques for hydrological simulations of climate change[END_REF][START_REF] Smitha | An improved bias correction method of daily rainfall data using a sliding window technique for climate change impact assessment[END_REF]. The size of a moving window can vary depending on the particular applications, e.g. [START_REF] Themeßl | Empirical-statistical downscaling and error correction of daily precipitation from regional climate models[END_REF] used a 61-day window to correct daily precipitation with the EQM method.

Long-term calibration

In another example, [START_REF] Switanek | Scaled distribution mapping: a bias correction method that preserves raw climate model projected changes[END_REF] calibrated climate data in smaller time blocks to better preserve the temporal evolution of climate change. They used 30-year periods with a 10-year sliding window to calibrate the ten middle years. For instance, the 2031-2040 period was calibrated using a model set up in the 2021-2050 (future) versus 1971-2000 (historical) periods. The users can adjust several parameters such as the length of the calibration period (30 years) and the length of the middle section (10 years). This process can be done using a sliding of, e.g. 10 years (to calibrate the 2041-2050 period). This sliding window technique can help the calibration better capture the raw modelled temporal evolution of climate change signal. However, it might not be suitable for all calibration methods, e.g. a direct method like QM calibrates the future data using the same transfer function. In this case, using the sliding window technique will not help gain more information.

Spatial configuration

Generally, the transfer function used to calibrate the data is "trained" at a pixel level, i.e. each modelled pixel is calibrated independently towards a corresponding observed pixel. These pixel-wise calibrated results, thus, are very close to the observations. However, there could be a problem of discontinuities when moving from one pixel to the next as -in this case -we have a distinct model for each pixel. A solution could be to consider a single model applicable for all pixels, but this could bring the calibrated data less close to the observations. We propose here a compromise that uses a spatial moving window: e.g. the 3 × 3 moving window shown in Fig. 6.10 uses one pixel of interest and eight surrounding pixels to train a local calibration model. Then, this localized transfer model is used to calibrate the central pixel.

Figure 6.11 compares the results of the pixel-wise versus the spatial moving window approaches for January temperature. The model data are calibrated and averaged over the future period (2076)(2077)(2078)(2079)(2080)(2081)(2082)(2083)(2084)(2085)(2086)(2087). The EqQM calibration method is chosen here, but other methods could be easily applied. At first glance, the two techniques seem to give similar results (Cal in the second column) as it is difficult to see a difference between the two calibrations. However, discrepancies can be noted when computing the (Cal -Mod) differences; in particular, the moving window approach provides smoother differences, meaning that the calibration process is smoother between the pixels. When zooming in one region (last column of Fig. 6.11), the differences in the calibration resolution becomes clear. By working at the pixel level, the impact of the observations is strong as the spatial pattern will be preserved in the calibrated data. This spatial pattern, especially at high resolution, is very much related to the land use pattern; and it is not possible to extrapolate this land-use pattern in the future. In the moving window approach, what is intended is to emphasize more the climate change signal from the model and give less emphasis on the observations. As often, this is a matter of choice, and the climate modellers need to decide based on what is important in its application. The size of the moving window is a parameter that can help them balance the two aspects.

Discussions and conclusions of the chapter

Calibration techniques become more and more important as they are essential to exploit climate change simulations. In particular, the climate impact studies measuring the socio-economic consequences of climate change require calibrated data in many fields (e.g. energy, agriculture, hydrology, and water management). A calibration technique often has two goals: (1) to make the climate simulations more realistic by bringing them closer to observations in a historical record; but at the same time, (2) to preserve the climate change signal that is provided by the model simulations. Each technique will make a particular compromise between these two aspects. The compromise will be dictated by the statistical properties that climate scientists want to obtain on the calibrated data. No method can be said to be better than the others; it is just a choice made on what is essential for the user. To facilitate this choice, a good understanding of each approach is necessary. Available calibration methods have become more complex, and it is more difficult for the user to choose and decide which one to use. Therefore, an up-todate synthetic review of available methods was necessary. This chapter tried to do so by presenting in a pedagogic way the available methods, emphasizing the main ideas, the mathematical formulas, the pertinent references, and their advantages and drawbacks.

It is not easy to judge if one method is better than another. For instance, it is possible to obtain a method that makes the climate simulations closer to the observations; however, this does not mean that the obtained transformation will be correct in the future, and no direct evaluation can be done. Thus, diagnostics need to check both the proximity of the simulations to the observations and the preservation of climate change information. Several diagnostics were introduced for that reason to measure the impact of using one technique over another, such as climate change signals in the mean, standard deviation, skewness, or trend.

For the three calibration approaches (i.e. direct, delta, and combined), experiments have been conducted for temperature and precipitation (two variables with a very different nature) over Europe. These experiments were considered using monthly data, but a daily scale would be necessary if the important statistical property was related to, for instance, the number of hot days or rain intensity. Our diagnostics can show the limitations and advantages of each approach. Overall, the combined approach seems to be a good candidate because it preserves the climate change signals well. Some aspects of the calibrated data can be improved. For example, the required calibrations could be constrained on a seasonal scale. Long-term tendencies from the model can also be better preserved using the moving window calibration centred on a decade. Furthermore, the balance of observed information versus model outputs can also be controlled; we proposed using the spatial resolution on which the transfer function of the calibration is set up for this purpose.

Numerous perspectives can be considered and discussed:

Calibration and downscaling -In the literature, the calibration of climate model simulations is often associated with downscaling (from a coarse resolution grid of the model to a higher resolution of the observations). Resolutions of 1 km can sometimes be attained. Simple regression models have been traditionally used to downscale model simulations in weather forecasting, and this has been naturally extended to climate simulations. Downscaling, in this context, means that the spatial patterns present in the historical record are projected to the future. This can be a good thing, but it can be dangerous too because high-resolution patterns might be related to, for instance, land use, which can be completely different in the future. Dissociating the downscaling and the calibration can therefore be beneficial. The spatial moving window approach that was proposed here might be a solution to balance the use in the calibration of this high-resolution information from the observation.

Calibration and model errors -Calibration can sometimes be used to correct important issues in the climate model. For example, a climate oscillation pattern in the model could not be well located compared to real observations. Using a calibration model to correct such deficiencies is too large a stretch of what can be asked to simple statistical models. It is illusory to ask to a simple regression model, at the pixels level, to fix a complex physical model. Therefore, calibration should not be used to pretend to correct such important climate features.

How to evaluate a calibration? -It is not easy to see if a calibration method is of good quality or not. For instance, it can perform well for the historical record but not be adequate for the future that cannot be evaluated. There is even the risk of getting too close to the observations to the detriment of what is done in the future (i.e. this issue is called "over-training" in machine learning theory). The solution proposed here is to use a large set of diagnostics showing the properties of the 133 original and the calibrated data and decide if these diagnostics are satisfactory for the application that has been considered.

Calibration of extreme values -

The calibration of extreme values is a challenging problem for several reasons. First, the definition of extremes is tricky: it can be based on quantile information (highest and lowest values at the 5 or 10 % level), a number of (rain/hot) days higher than a threshold, a frequency of occurrence of a particular event, etc. Each extreme definition would require a particular diagnostic tool, a particular transfer function for the calibration. For particular aspects of extreme cases, hand-tailored models need to be designed. Second, a statistical method such as a calibration is mainly driven by the most common events in the database used to set it up, not by the rarest extreme events. Therefore, the applicability of the calibration model to extreme values can be questionable and erroneous. Third, future extreme values in climate simulations can extend beyond the observed ranges in the historical record used to set up the calibration. This out-of-range will jeopardise the applicability of the calibration model for future scenarios. Delta methods utilise relative changes and can overpass this problem, but direct and combined approaches can suffer from such issues. In this case, some form of extrapolation, e.g. proposed by [START_REF] Wood | Hydrologic implications of dynamical and statistical approaches to downscaling climate model outputs[END_REF][START_REF] Themeßl | Empirical-statistical downscaling and error correction of regional climate models and its impact on the climate change signal[END_REF], could help deal with these outside values. Finally, extremes are even more vivid for climate simulations on a daily scale, as the extreme definition is even more complex and comparing observations and model outputs will be very challenging.

Multivariate calibration -

The calibration methods proposed in this chapter deal with marginal distributions, meaning that the adjustment was made independently in each considered variable (e.g. temperature or precipitation), ignoring the link between/among variales and influence of all other variables (or locations) [START_REF] Maraun | Nonstationarities of regional climate model biases in european seasonal mean temperature and precipitation sums[END_REF]. This approach risks breaking the link between the climate variables and introducing physical inconsistencies. A good measure of the dependency among climatic variables could help reduce this problem [START_REF] Dosio | Bias correction of the ensembles high-resolution climate change projections for use by impact models: Evaluation on the present climate[END_REF]. Several studies recently suggested using a multivariate method dealing with the joint distributions of the climate variables [START_REF] Maraun | Bias correcting climate change simulations -a critical review[END_REF]. Such methods aim at adjusting the joint distributions without breaking the dependency structure between the variables [START_REF] Vrac | Multivariate-intervariable, spatial, and temporal-bias correction[END_REF][START_REF] Cannon | Multivariate bias correction of climate model output: Matching marginal distributions and intervariable dependence structure[END_REF][START_REF] Vrac | Multivariate bias adjustment of high-dimensional climate simulations: the rank resampling for distributions and dependences (r 2 d 2 ) bias correction[END_REF]. However, most multivariate methods are bivariate as they conserve the dependency between two variables only. More variables could theoretically be considered when adjusting joint distributions, but ultimately, all the climate variables should be included. In this case, the statistical model would become a full climate model, and this ambitious task is obviously not possible. This multivariate aspect is a fundamental limitation of the calibration principle: between a practical but very simplistic transformation and a more complex statistical model that is unattainable. This chapter finally focuses on the climate change impacts on coffee production. As in the climate change context, many factors, such as management practices and agriculture policies, might change dramatically. Therefore, estimating the yield in the future might not be practical. Here, we are less ambitious compared to the purposes of two previous chapters (chapters 4 and 5). This chapter aims to predict how suitability for coffee production should change in the future with an assumption of equilibrium, meaning that coffee plants do not adapt to new conditions under climate change and no adaptation strategies are implemented (e.g. new varieties, management of microclimate with shade trees, etc.). As a less ambitious task, the suitability model requires fewer data (i.e. coffee data). Consequently, we can perform the test in more coffee regions, even when there are only about a few years of data. Examples of the two largest coffee regions are presented in this chapter: the Arabica-growing areas in Brazil and the Robusta-growing areas in Vietnam. Climate data, calibrated by a combined approach (i.e. EqQM presented in Sect. 6.2.2.3), will be used to derive the model inputs. Although farm-based data were used in previous studies (e.g. in Bunn et al. [2015b,a]), these data are actually quite complex and need to be systematic. We will propose a more direct approach to determining coffee occurrence data derived from regional/national statistics on coffee-growing areas.

Artificial intelligence (AI) in calibration -

This chapter is organized as follows. Section 7.1 introduces the needed materials such as the study area, coffee data, potential predictors, and climate data. Then, the methods used in this chapter are described in Sect. 7.2. Section 7.3 presents the main results of the model validation and of the suitability model (e.g. the current and future suitabilities, the sensitivity to emissions scenarios, and the distribution of climate change impacts). Discussions and conclusions are in Sect. 7.4.

This thesis chapter is the subject of the manuscript submitted to Earth's Future entitled "Climate change impacts on Robusta coffee production in Vietnam." [START_REF] Dinh | Climate change impacts on robusta coffee production in vietnam[END_REF].

Materials

Study area and coffee data

Two areas (shown in Fig. 7.1) will be considered in this chapter.

• (1) Arabica area -This study area is located along the longitude 57.5 • W to 37 • W and latitude 12 • S to 28 • S, which covers the Brazilian high-intensity coffee production area. In this area, most of the produced coffee is Arabica (more than 90 %); we will thus consider it as an Arabica example.

A map of the study area is shown in Figs. 7.1a1-c1. In detail, Fig. 7.1a1 shows the South America domain covered by the regional climate models (i.e. CORDEX-SAM22, see in Sect. 7.1.3). Then, Figs. 7.1b1 and c1 provide a close-up of our Arabica study area with the information on the topographical elevation and soil type, respectively. Digital elevation data are downloaded from the USGS EROS Center, and soil type data are from the FAO website, as presented in Sect. 2.1.1.4.

• (2) Robusta area -The second study area is located along the longitude 104.22

Methods

Learning database identification

To develop a suitability model, we need a database of samples indicating the presence (coffee samples) or absence (background samples) classes. Such a database can therefore be used to calibrate a statistical suitability model. The presence data are based on observations (i.e. data from GSO); the absence data are generated and not observed. The detailed methodology to identify these two data is presented as follows.

The current coffee samples are built from the available coffee areas with respect to the climate-gridded cell. Here, we neglected coffee-growing districts with less than 500 ha as they are relatively small compared to the size of one climate cell. For municipalities/districts with coffee-growing areas of 500 ha, one corresponding gridded cell will be selected as a coffee sample (or coffee cell): (1) if the shape of the municipality/district is smaller than the cell or it contains only one cell, we chose that gridded cell; (2) if the shape is overlapped by several cells (e.g. often two cells), the cell with higher overlap surface was selected; (3) if the shape contains several cells which was seldom the case, the coffee cell was randomly chosen among these cells. For municipalities/districts with coffee-growing areas higher than 500 ha, several corresponding cells were considered and classed as current coffee samples. For each municipality/district, the number of cells was defined by a scaling factor S f , which is computed by the ratio of the coffee area over 500 ha. For instance, let's consider a municipality/district with 5000 ha of coffee. The corresponding scaling factor S f = 5000 500 = 10, meaning that there will 141 be ten coffee cells associated to this municipality/district. We identified these ten cells by replicating the corresponding gridded cell(s).

The background samples are randomly selected from a 4.4 • buffer around present regions. The ratio of background samples to current coffee data is set to 1:1, as recommended in [START_REF] Barbet-Massin | Selecting pseudo-absences for species distribution models: how, where and how many?[END_REF]. Too few (or too many) background samples can lead to overfitting (or underfitting) as the model is biased toward the coffee samples (or the background samples).

The learning database was used to explore the corresponding elevation and soil information of the coffee and background samples. This information was used to exclude points that are unsuitable for coffee for both current climatic conditions and future climate projections.

For Arabica coffee, Fig. 7.2a1 shows the map of the learning database -coffee (red points) and background (blue points) samples -together with the coffee area data over the study area. We also plot the normalised histogram of the elevation (in m) and the soil type corresponding to these coffee and background points in Figs. 7.2b1 and c1, respectively. The learning database suggests that the Arabica coffee requires an elevation higher than 200 m, while most coffee samples range from 600 to 1200 m. Also, Arabica grows well in clay-loam, clay, sandy-clay-loam, or loam soils. Among the selected coffee samples, a very small percentage of coffee points is grown in loamy-sand, sandy-clay, or sandy-loam soils. No coffee is grown in sand, heavy-clay soils, as well as water. Similarly, Fig. 7.2a2 to c2 show the learning database results for Robusta coffee. Figure 7.2b2 suggests that Robusta coffee requires an elevation higher than 100 m, and the optimal range is from 500 to 800 m. Concerning the soil type, Robusta grows well in clay, clayloam, or sandy-loam soils. In contrast to Arabica, a very small percentage of coffee points is grown in sandy-clay-loam soil, and no coffee is grown in loam soils. Other soil types are not found in this Southeast Asia region. In the climate change impact assessments, we excluded all cells below 200 m and/or sand, heavy-clay soils for Arabica coffee, and the cells below 100 m and/or loam soil for Robusta coffee.

Suitability model

Model selection

Our suitability model relies on neural networks (NN) (as described in Sect. 3.1.3). As in a classification task, NN trains the generic feedforward neural networks to map each input vector into its corresponding target vector. The target is based on the learning database introduced previously. The inputs are chosen from 12 potential predictors (Sect. 7.1.2). A forward selection method (Sect. 3.5.1) is used to obtain the hierarchy of the explanatory variables. However, the final model only

Impacts

We first trained and tested the model on the learning database for the current climate (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019) to obtain optimal model parameter values. Then, the model is applied to all cells in the study area for different periods: current (i.e. 2009-2019) and future (i.e. 2038-2048, 2049-2059, and 2060-2070 periods). The results are visualized as maps with continuous scores, which are normalised from 0 (not suitable) to 1 (suitable). A threshold is used to determine if a cell is suitable or not for coffee. The threshold is based on the probability density functions of these two classes in the actual learning database (Fig. 7.3) and the coffee area. As shown in Fig. 7.3a and b, a threshold of 0.5 helps to distinguish very well between the coffee and background samples for both Arabica and Robusta coffee cases. In addition, with this threshold, our model was able to identify about 1 000 × 10 4 ha of suitable areas in the Arabica area and approximately 70 × 10 4 ha of suitable areas in the Central Highlands of Vietnam. These values are comparable to the actual coffeegrowing area in these two regions. The area is computed by summing the entire land area within suitable cells (i.e. a pixel cell of 0.22 • × 0.22 • ). For suitability assessments of future climate projections, we first applied the suitability model for each of the calibrated RCMs shown in Tab. 7.2 (one for Arabica coffee and eight for Robusta coffee). Then, we compared the suitability changes under RCP2.6 and RCP8.5 scenarios using GERICS_NCC_0 model for Arabica coffee and six RCMs (i.e. ICTP_Had, GERICS_Had, ICTP_NCC, GER-ICS_NCC, ICTP_MPI, and GERICS_MPI) for Robusta coffee. The SMHI_CNRM and SMHI_Had models are excluded as they do not provide simulations for the RCP2.6 scenario. We finally investigated the area changes respective to elevation.

Calibration of climate simulations

As described in chapter 6, the combined approach (Sect. 6.2.2.3) seems to be a good candidate as it preserves the climate change signals well. Here, we used the equidistant quantile mapping (EqQM) [START_REF] Li | Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change ar4 models using equidistant quantile matching[END_REF][START_REF] Pierce | Improved bias correction techniques for hydrological simulations of climate change[END_REF][START_REF] Wang | Equiratio cumulative distribution function matching as an improvement to the equidistant approach in bias correction of precipitation[END_REF]] to calibrate the RCMs (i.e. CORDEX-SAM and CORDEX-SEA data). The calibration was done using three datasets, including:

• the historical observations X O,h (i.e. the integrated ERA5-Land from 1981 to 2005),

• the corresponding simulations X M,h on the historical record (i.e. CORDEX-SAM22 and CORDEX-SEA22 from 1981 to 2005),

• and the simulations X M,f for the future (i.e. RCMs for three periods: 2031-2055, 2042-2066, and 2053-2077).

In addition, the calibration is done for several time blocks to preserve the temporal evolution of climate change better (Sect. 6.4.1). Here, we used 25-year periods with an 11-year sliding window to calibrate the 11 middle years. For instance, the 2038-2048 period was calibrated using a model set up in the 2031-2055 (future) versus 1981-2005 (historical) periods. Next, we calibrated 2049-2059 data using the future 2042-2066 period and the same historical period. Finally, we used the model set up in the 2053-2077 and 1981-2005 periods to obtain the 2060-2070 calibration.

Results

Model validation and variable contribution

The suitability models, which are trained and tested over the learning database, show a good predictive performance for the current period (2009-2019), both Arabica and Robusta coffee. For the Arabica case, the model precision is 92.3 %, meaning that less than 8 % of background samples are misclassified as coffee regions. We also obtain a high model recall of 89.1 %, implying that the model can identify the actual coffee cells well on the historical record based on actual observations. These results are even better for Robusta, with a precision of 95.4 % and a very high model recall of 96.9 %.

Among 12 potential predictors, the mean maximum temperature during the growing season (i.e. Tmax39 shown in Tab. 7.1) was found to be the most important for Arabica. This was followed by the precipitation seasonality (Bio15) and the maximum number of consecutive dry months (NDM). Precipitations during flowering (P13) and during the growing season (P39) ranked as least important for our Arabica case. For Robusta, the most important variables are (1) the mean minimum temperature during the harvest and growing season (i.e. Tmin1012 and Tmin39); (2) the precipitation variables are also crucial for Robusta coffee, including precipitation during the late growing season (P79) and flowering (P13). In contrast to Arabica, the least important variables are the mean maximum temperature during the growing season (Tmax39) and the maximum temperature in the warmest months (Bio5).

Climate suitability 7.3.2.1 Current and future suitabilities

For Arabica, the (trained and tested) suitability model is applied over all study cells for the current period, i.e. 2009-2019, as shown in Fig. 7.4a. The highly suitable areas are in the Southeast region of Brazil (e.g. Minas Gerais state), which is in agreement with the current coffee planted area shown in Fig. 7.2a1. We then 2038-2048, 2049-2059, and 2060-2070, respectively. As in the peak scenario, the total suitable area gradually reduces and reaches a very low number in the 2054s. After that, this number starts to increase in the 2065s.

Similar to the Arabica study, the suitability model is first applied for the 2009-2019 period, and then future suitabilities are investigated for different periods. However, eight RCMs (i.e. models No. 2 to 9 presented in Tab. 7.2) are considered instead, and we will consider the high greenhouse gas concentration scenario (RCP8.5). As shown in Fig. 7.5a, the largest current suitable areas are in the mountainous locations, especially the Central Highlands of Vietnam. Again, these results are similar to the current coffee planted area shown in Fig. 7.2a2. For . 7.5b2, c2, and d2. Compared to the current suitability, highly suitable areas will decrease significantly in the future. These areas become smaller and smaller in time. For example, the suitable area is about 81 × 10 4 ha for the current period (Fig. 7.5a); however, it will be reduced by half in the 2065s, i.e. ≈ 42 × 10 4 ha (Fig. 7.5d1). In addition, the eight simulations give similar predictions by showing small standard deviation values, i.e. 15 % on average, for all three projected periods.

Sensitivity to climate scenarios

We now study the suitability model's sensitivity to different climate scenarios. For Arabica, results are induced from only one regional model (i.e. GERIS_NCC_0). For Robusta, on the other hand, we will investigate the ensemble results from six RCMs (i.e. model No. 4 to 9 in Tab. 7.2) since they provide the simulations for both RCP2.6 and RCP8.5 scenarios. First, let us look at the suitability changes (the difference between future and current suitability) for different future periods under the RCP2.6 and RCP8.5 scenarios shown in Fig. 7.6, for Arabica. As expected, climate change impacts are less pronounced in the low CO2 emissions scenario (a1 to a3) than in the high CO2 emissions scenario (b1 to b3). For the 2038-2048 period, for instance, the number of magenta cells, which signify the negative change, is much less and lighter in (a1) than (b1). Similar behaviours are obtained for two other considered periods, as shown in (a2) versus (b2) and (a3) versus (b3). In addition, both scenarios show that most of the current growing Arabica regions suffer negative impacts or could become unsuitable, i.e. shown in many magenta cells in Fig. 7.6. Nevertheless, many areas in the southern margin of Brazil (e.g. Santa Catarina, Parana, Sao Paulo, and Rio de Janeiro) could become more suitable than present for arabia. For the RCP8.5 scenario, the negative impacts become stronger toward the end of the century: more regions suffer from negative impacts in the 2060-2070 period (b3) than in the 2038-2048 period (b1). Figure 7.7 shows the suitability changes for Robusta over three future periods under the RCP2.6 and RCP8.5 scenarios. Similar to Arabica, the Robusta production regions may see severe losses of suitability, with possible positive changes in the northern Central Highlands. This area covers high mountains ranging from 875 to 1200 m. The negative impacts are observed to be more significant under the RCP8.5 scenario than those in the RCP2.6 scenario. With this high-emissions scenario, more and more regions suffer from negative changes toward the end of the century, e.g. many dark magenta cells in Fig. 7.7b3. Figures 7.8 and 7.9 show the distribution of the suitable regions by elevation for three future periods (i.e. 2038-2048, 2049-2059, and 2060-2070), for Arabica and Robusta, respectively. Here, we compute the suitable areas by summing the entire land area within suitable cells across 100 m elevation classes. For Arabica, results from the GERICS_NCC_0 model suggest the loss of suitability, mostly in low elevations, from 500 to 1 200 m for RCP2.6 and below 1 200 m for RCP8.5. For the RCP2.6 scenario, the largest loss of total suitable area is in the 2049-2059 period, in which the total suitable area decreases by 18 %. On the other hand, for the RCP8.5 period, the total suitable area decreases gradually in time, and losses are projected to be 34 % by the mid-2060s. More results can be interpreted for Robusta, including the ensemble mean and uncertainties (shaded areas shown in Fig. 7.9) provided by several RCMs. In general, suitable areas under the RCP2.6 scenario (green lines and shaded areas in Fig. 7.9) are somewhat comparable to the currently suitable areas (black lines in Fig. 7.9). The average losses are about 20 % over three projected periods. The most considerable loss of suitable area could be up to 36 % for the 2060-2070 period. For the high impact scenario RCP8.5, the suitable areas decrease significantly for three projected periods, with the losses ranging from 39 % up to 83 %. The total suitable areas do not change very much from 2038-2048 to 2049-2059; however, we observe a substantial decrease in the mid-2060s, i.e. about 15 % compared to the two previous projected periods. The suitability is shown to be very sensitive to the elevation. In both scenarios, significant losses are found in low-elevation areas (i.e. below 800 m). In contrast, the total suitable areas do not change much in higher elevations, especially above 850 m.

Discussions and conclusions of the chapter

This chapter assessed the climate change impacts on coffee production in the two largest coffee-growing regions (i.e. Arabica coffee in Brazil and Robusta coffee in Vietnam). The suitability models indicate that projected climate change scenarios will negatively impact the suitability for growing coffee in the two studied areas, even in the low-impact scenario (RCP2.6). The degree of suitability change depends on the emissions scenarios and periods. Overall, the losses are higher for Robusta than for Arabica. As expected, suitability decreases towards the end of the century, particularly after 2060. Results from one climate model suggest losses noticeably in suitability for Arabica in Brazil, especially for the RCP8.5 scenario. Using an ensemble of models will probably help capture more interesting information. For Robusta, significant losses are found in the Robusta-growing area in Vietnam, and this is in line with a previous study of Bunn et al. [2015b]. These losses could probably be due to the increased temperature in the climate change context since the temperature (particularly minimum temperature) is a major determinant of coffee's distribution of climate suitability.

Both species are sensitive to elevations: the losses in suitability are mostly at low elevations; conversely, the suitability increases at higher elevations. At higher elevations, coffee can benefit from less increase in temperature within the climate change context. These results are consistent with previous studies focusing on local and global coffee-growing regions [START_REF] Schroth | Towards a climate change adaptation strategy for coffee communities and ecosystems in the sierra madre de chiapas, mexico[END_REF]Bunn et al., 2015b]. Mitigation strategy by shifting to higher elevations will likely be a global trend, even though further tests are needed to identify the magnitude of this impact. In addition, fewer suitable areas for coffee (or crop in general) are available at higher elevations because these sites are often quite steep and have low permanent soil cover. This limitation will make mitigation strategies more challenging.

We develop here a new technique to define the coffee sample (i.e. present occurrence data) based on the coffee area data. This dataset is direct and easier to obtain compared to the classical farm-based data. It thus allows us to work globally where the coffee information is available, even for a few years of data. Our identification technique here is very simple and assumes that the current coffee distribution is an appropriate sample for suitable growing conditions, but these conditions are not only affected by climate. In future works, other factors, such as the land-use classes to better identify the coffee presence locations, thereby improving the quality of the suitability models.

Due to the time limitations of the thesis, this chapter focuses only on two main areas, and only one model is included in the Arabica study. Pursuing this work, we can apply the suitability model to other coffee regions (e.g. the ones presented in Sect. 2.1.1.3) since fewer coffee data (e.g. coffee areas) are required in this application. It is also interesting to have an ensemble of climate models across several emissions scenarios instead of a single model. Also, methodologies developed in this chapter can be easily applied to other crops and regions.

Conclusion

In this thesis, we develop and test several statistical approaches to understand the impact of weather and climate on crop yields, especially for coffee. The first approach uses a simple linear regression model with regularization techniques and leave-one-out cross-validation. The LOO method is commonly used; however, there are some limitations, and it can also be misused in several cases, especially when having a small number of samples. Therefore, we develop the LTO method to choose a suitable crop model and assess its true generalization quality.

Two approaches are applied to estimate the Vietnamese Robusta coffee yield. Results suggest that temperature and precipitation information can account for up to 36 % of yield anomalies, which is in agreement with previous studies. This sensitivity of yield anomalies to weather can vary substantially among Vietnamese provinces and districts due to the large spatial weather conditions of the studied regions. Also, the correlation analysis suggests that Arabica in Brazil is less weather-sensitive than Robusta in Vietnam. Results show that Vietnamese Robusta coffee is most sensitive to two key moments: a prolonged rainy season of the previous year favouring vegetative growth, thereby increasing the number of fruiting nodes and the potential yield, while low rainfall during bean formation decreases yield. Overall, the model can forecast Robusta coffee yield anomalies from three to six months before the harvest.

In addition, to demonstrate that the LTO method is general and easily extended to other regions and crops, we also perform the test on French grain maize -a very mature culture. Estimates of grain maize yield at the end of the season show that weather can explain more than 40 % of the variability in yield anomaly. The LTO method is also useful in seasonal yield forecasting: mid-season grain maize yield forecasts are possible from June, and weather variables in July are among the most influential predictors.

Climate simulations are our unique sources of information for climate change; however, they are generally biased and cannot be used directly in climate change impact studies. Numerous methods to correct/calibrate climate model simulations have recently been developed. This thesis provides an up-to-date synthetic review of the available calibration methods. These methods can be divided into three approaches: direct, delta, and combined. Each method will make a particular compromise between two aspects: (1) to make the climate simulations more realistic by bringing them closer to observations in a historical record; (2) to preserve the climate change signal that is provided in the model simulations. No method can be said to be better than others; it is just a choice of what is essential for the user. We choose here to use the combined approach based on thorough diagnostics. We propose new and simple diagnostic tools to compare different methods, and techniques to improve spatio-temporal improvements configurations of the calibration are also discussed.

For the long-term assessments, we focus rather on estimating the change in the coffee suitability area instead of quantifying the change in future yield, as many factors can change in the future context. We use the calibrated climate simulations and suitability models that include both current coffee data and static information (soil and elevation) to estimate the changing area over two important coffee-producing areas (i.e. the Arabica-producing area in Brazil and the Robustaproducing area in Vietnam). Results show that projected climate change scenarios will negatively impact the coffee suitability areas. The degree of suitability change depends on the emissions scenarios and periods: the suitable Arabica areas in Brazil could decrease by about 26 % by the mid-century in the high-emissions scenario, and this number is surprisingly high for Vietnamese Robusta (≈ 60 %). Significant losses are found in elevations below 1 200 m for Arabica and 800 m for Robusta, suggesting a possible adaption strategy by shifting to higher locations.

Discussions and perspectives

There are numerous perspectives for pursuing this thesis. For instance, concerning the methodology, a (crop) model is generally structured into three parts: the inputs, the model, and the outputs. Each section can be improved to obtain better crop models. First, our statistical models mainly considered the weather information, even though our suitability models included elevation and soil types as time-invariant conditions. In future works, we could also include other predictors such as agro-climatic indicators characterising plant-climate interactions for global agriculture or other relevant predictors (e.g. (very) heavy precipitation days or frost days, for coffee). The addition of input predictors should not come at the cost of over-training, as more input information is possible only when more data (i.e. yields) is available for the calibration of the model. Then, mixed-effects models -which use samples in several regions -can be considered to compensate for the lack of historical data. Similar to the idea of the mixed-effects model, a multi-crop model can also be feasible. By having groups of different crops, this model can take into account the specificities of each group (i.e. crop type) and the variability between groups. Such a model could benefit from a more extensive database, potentially exploit the crops' dependency, and facilitate risk management on a farm, national, or even global scale. For instance, model results can suggest farmers adopt the practice of multiple-cropping instead of mono-cropping. This system can, for example, reduce the risk of crop failure. Other models, such as machine learning and deep learning (e.g. (convolutional) neural networks), can also be used, but this would require more samples. In addition, a combining approach of statistical and process-based models can also be of interest for regional study as it inherits the advantages of the two model types. This combined approach has been shown to be remarkably better than either model [START_REF] Roberts | Comparing and combining process-based crop models and statistical models with some implications for climate change[END_REF][START_REF] Gornott | Improving crop modeling approaches for supporting farmers to cope with weather risks[END_REF]. Here, for instance, we can use a process-based model as a regularization tool to identify the right dependency structure of a statistical model.

Throughout the thesis, our statistical models are applied in three modes: yield forecasts during the season (i.e. seasonal forecasts), end-of-season yield estimations, and long-term predictions. These applications are very important and helpful in management, monitoring, or defining adaptation strategies. For example, in the current crisis, having quick (and reasonable) estimations of crop yields is essential for many countries, especially those that highly rely on Ukraine agriculture exports, to make decisions regarding stock management, imports and exports. In future works, we can extend our statistical models and their applications to other cultures and regions in the world, e.g. other French crops (winter/spring soft/durum wheat, barley, oats, etc.), soft/durum wheat over 27 EU countries, or coffee in Peru. Real-time data from satellites can be involved in these tasks, but these data are only valid for historical and current times. In addition, it would be interesting to test the model with finer temporal resolution (e.g. weekly or daily) to provide more information during the growing season. Like coffee, climate change impact assessments can be performed for French crops. Other GCMs or RCMs for other scenarios can also be considered for coffee and other cultures. Finally, the purpose of the crop models can be fulfilled if their outputs become more accessible to different users (farmers, local authorities or ministries). This could be done by developing an online platform with the help of experts in the field to make quick decisions during the growing season and before/after the harvest.

La météo et le climat ont un impact important sur les rendements agricoles. De nombreuses études basées sur différentes approches ont été réalisées pour mesurer cet impact. Cette thèse se concentre sur l'utilisation de modèles statistiques pour mesurer la sensibilité des cultures aux conditions météorologiques sur la base des enregistrements historiques. Lors du développement et de l'utilisation d'un modèle statistique, une difficulté critique survient lorsque les données sont rares, ce qui est souvent le cas pour la modélisation des cultures. Il y a un risque élevé de sur-apprentissage si le modèle n'est pas développé avec certaine précautions. Ainsi, la validation et le choix du modèle sont deux préoccupations majeures de cette thèse.

Deux approches statistiques sont développées. La première utilise la régression linéaire avec régularisation et validation croisée (c.-à.-d. leaveone-out ou LOO), appliquée au café robusta dans la principale région productrice de café du Vietnam (c.-à.-d. les Montagnes Centrales). Le café est une culture rémunératrice, sensible aux intempéries, et qui a une phénologie très complexe en raison de sa nature pérenne. Les résultats suggèrent que les informations sur les précipitations et la température peuvent être utilisées pour prévoir l'anomalie de rendement avec une anticipation de 3 à 6 mois selon la région. Les estimations du rendement du robusta à la fin de la saison montrent que les conditions météorologiques expliquent jusqu'à 36 % des anomalies de rendement historiques. Cette première approche de validation par LOO est largement utilisée dans la littérature ; cependant, elle peut être mal utilisé pour de nombreuses raisons : elle est technique, mal interprétée et nécessite de l'expérience. Une alternative, l'approche "leave-two-out nested cross-validation" (ou LTO), est proposée pour choisir le modèle approprié, évaluer sa véritable capacité de généralisation et choisir la complexité du modèle optimale. Cette méthode est sophistiquée mais simple. Nous démontrons son applicabilité pour le café robusta au Vietnam et le maïs en France. Dans les deux cas, un modèle plus simple avec moins de prédicteurs potentiels et d'entrées est plus approprié. Utiliser uniquement la méthode LOO, sans aucune régularisation, peut être très trompeur car cela encourage à choisir un modèle qui sur-apprend les données de manière indirecte. L'approche LTO est également utile dans les applications de prévision saisonnière. Les estimations de rendement du maïs en fin de saison suggèrent que les conditions météorologiques peuvent expliquer plus de 40 % de la variabilité de l'anomalies de rendement en France.

Les impacts du changement climatique sur la production de café au Brésil et au Vietnam sont également étudiés à l'aide de simulations climatiques et de modèles d'impact (appelées "suitability models"). Les données climatiques sont cependant biaisées par rapport au climat réel. Par conséquent, de nombreuses méthodes de "correction de biais" (appelées ici "calibration") ont été introduites pour corriger ces biais. Une présentation critique et detaillée de ces calibrations dans la littérature est fournie pour mieux comprendre les hypothèses, les propriétés et les objectifs d'application de chaque méthode. Les simulations climatiques sont ensuite calibrées par une méthode basée sur les quantiles avant d'être utilisées sur nos modèles d'impact. Ces modèles sont développés sur la base des données de recensement des zones caféières, et les variables climatiques potentielles sont basées sur un examen des études précédentes utilisant des modèles d'impact pour le café et des recommandations d'experts. Les résultats montrent que les zones propices à l'arabica au Brésil pourraient diminuer d'environ 26 % d'ici le milieu du siècle dans le scénario à fortes émissions, les régions compatibles avec la culture du robusta vietnamien pourraient quant à elle diminué d'environ 60 %. Les impacts sont significatifs à basse altitude pour les deux types de café, suggérant des déplacements potentiels de la production vers des endroits plus élevés.

Les modèles statistiques que nous construisons, notamment la technique LTO, peuvent contribuer au développement de la modélisation des cultures. Ils peuvent être appliqués à une culture pérenne complexe comme le café ou sur des cultures annuelles plus industrialisées comme le maïs. Ils peuvent être utilisés dans les prévisions saisonnières ou les estimations en fin de saison, qui sont utiles pour la gestion et le suivi des cultures. Estimer l'adéquation future des cultures permet d'anticiper les conséquences du changement climatique sur le système agricole et de définir des stratégies d'adaptation ou d'atténuation. Les méthodologies utilisées dans cette thèse peuvent être facilement généralisées à d'autres cultures et régions du monde.

Weather and climate strongly impact crop yields. Many studies based on different techniques have been done to measure this impact. This thesis focuses on statistical models to measure the sensitivity of crops to weather conditions based on historical records. When using a statistical model, a critical difficulty arises when data is scarce, which is often the case with statistical crop modelling. There is a high risk of overfitting if the model development is not done carefully. Thus, careful validation and selection of statistical models are major concerns of this thesis.

Two statistical approaches are developed. The first one uses linear regression with regularization and leave-one-out cross-validation (or LOO), applied to Robusta coffee in the main coffee-producing area of Vietnam (i.e. the Central Highlands). Coffee is a valuable commodity crop, sensitive to weather, and has a very complex phenology due to its perennial nature. Results suggest that precipitation and temperature information can be used to forecast the yield anomaly with 3-6 months' anticipation depending on the location. Estimates of Robusta yield at the end of the season show that weather explains up to 36 % of historical yield anomalies. The first approach using LOO is widely used in the literature; however, it can be misused for many reasons: it is technical, misinterpreted, and requires experience. As an alternative, the "leave-two-out nested cross-validation" (or LTO) approach, is proposed to choose the suitable model and assess its true generalization ability. This method is sophisticated but straightforward; its benefits are demonstrated for Robusta coffee in Vietnam and grain maize in France. In both cases, a simpler model with fewer potential predictors and inputs is more appropriate. Using only the LOO method, without any regularization, can be highly misleading as it encourages choosing a model that overfits the data in an indirect way. The LTO approach is also useful in seasonal forecasting applications. The end-of-season grain maize yield estimates suggest that weather can account for more than 40 % of the variability in yield anomaly. Climate change's impacts on coffee production in Brazil and Vietnam are also studied using climate simulations and suitability models. Climate data are, however, biased compared to the real-world climate. Therefore, many "bias correction" methods (called here instead "calibration") have been introduced to correct these biases. An up-to-date review of the available methods is provided to better understand each method's assumptions, properties, and applicative purposes. The climate simulations are then calibrated by a quantile-based method before being used in the suitability models. The suitability models are developed based on census data of coffee areas, and potential climate variables are based on a review of previous studies using impact models for coffee and expert recommendations. Results show that suitable Arabica areas in Brazil could decrease by about 26 % by the mid-century in the high-emissions scenario, while the decrease is surprisingly high for Vietnamese Robusta coffee (≈ 60 %). Impacts are significant at low elevations for both coffee types, suggesting potential shifts in production to higher locations.

The used statistical approaches, especially the LTO technique, can contribute to the development of crop modelling. They can be applied to a complex perennial crop like coffee or more industrialized annual crops like grain maize. They can be used in seasonal forecasts or end-ofseason estimations, which are helpful in crop management and monitoring. Estimating the future crop suitability helps to anticipate the consequences of climate change on the agricultural system and to define adaptation or mitigation strategies. Methodologies used in this thesis can be easily generalized to other cultures and regions worldwide.

  les objectifs d'application de chaque méthode. Les simulations climatiques sont ensuite calibrées par une méthode basée sur les quantiles avant d'être utilisées sur nos modèles d'impact. Ces modèles sont développés sur la base des données de recensement des zones caféières, et les variables climatiques potentielles sont basées sur un examen des études précédentes utilisant des modèles d'impact pour le café et des recommandations d'experts. Les résultats montrent que les zones propices à l'arabica au Brésil pourraient diminuer d'environ 26 % d'ici le milieu du siècle dans le scénario à fortes émissions, les régions compatibles avec la culture du robusta vietnamien pourraient quant à elle diminué d'environ 60 %. Les impacts sont significatifs à basse altitude pour les deux types de café, suggérant des déplacements potentiels de la production vers des endroits plus élevés. Les modèles statistiques que nous construisons, notamment la technique LTO, peuvent contribuer au développement de la modélisation des cultures. Ils peuvent être appliqués à une culture pérenne complexe comme le café ou sur des cultures annuelles plus industrialisées comme le maïs. Ils peuvent être utilisés dans les prévisions saisonnières ou les estimations en fin de saison, qui sont utiles pour la gestion et le suivi des cultures. Estimer l'adéquation future des cultures permet d'anticiper les conséquences du changement climatique sur le système agricole et de définir des stratégies d'adaptation ou d'atténuation. Les méthodologies utilisées dans cette thèse peuvent être facilement généralisées à d'autres cultures et régions du monde.

Figure

  FIGURE 2.2: Coffee area (in ha) in 2017 for: (a) first and (b) second administrative levels.

FIGURE 2 . 4 :

 24 FIGURE 2.4: Coffee yield time series (in kg•ha -1 ) of the first and second administrative levels in Brazil and Vietnam: (a1) five major coffee-producing states in Brazil and (b1) five major coffee-producing provinces in Vietnam; (a2) 100 municipalities in Minas Gerais (Brazil) and (b2) several districts in Dak Lak (Vietnam).

  FIGURE 2.5: Elevation data over (a) Brazil and (b) Vietnam regions.

FIGURE 2 . 6 :

 26 FIGURE 2.6: Soil data over (a) Brazil and (b) Vietnam regions.

FIGURE 2 . 7 :

 27 FIGURE 2.7: Map showing 96 French departments and their corresponding names.

FIGURE 2 . 9 :

 29 FIGURE 2.9: Total wheat yield data (in t•ha -1 ): (a) yield in 2015 over 27 EU countries, and (b) yield time series of 16 German regions at NUTS 1 level.

FIGURE 2 .

 2 FIGURE 2.10: (a) Global annual total precipitation (in 10 3 mm), (b) global annual average 2 m temperature (in • C), in 2000.

FIGURE 2 .

 2 FIGURE 2.11: (a) Annual total precipitation (in mm) and (b) average annual temperature (in • C) over all the districts in the Central Highlands (Vietnam), averaged from 2010 to 2018.

  FIGURE 2.13: Description of the two considered CORDEX domains, adopted from WCRP [2015]: (a) South America (SAM) and (b) Southeast Asia (SEA). These regions are defined by several parameters using a rotated pole coordinate: (1) coordinates of the rotated pole in rotated coordinates: RotPole (Longitude; Latitude); (2) coordinates of the top left corner (TLC) in rotated coordinates: TLC (Longitude; Latitude); (3) number of grid point in the East-West direction: Nx;(4) number of point in the North-South direction: Ny.

  FIGURE 3.1: A comparison between (a) a biological neuron and (b) a NN neuron. In the NN neuron, synapses are modelled by the set of inputs x i with i = 1, 2, • • • , n. The cell body is modelled by the biological counterpart functionality, that is in collecting together weighted inputs and filter them throughout an activation function f . The figure is adopted from Roffo [2017].

FIGURE 3

 3 FIGURE 3.2: Structure of a NN model with n input inputs, one hidden layer with n neuron neurons, and one output layer with a singer output.

  FIGURE 3.3: (a) Graphical illustration of the bias and variance. (b) Schematic illustration of the bias-variance trade-off. Model bias decreases with increasing model complexity while model variance increases with increasing model complexity (adopted from Hastie et al. [2009a]).

FIGURE 3 . 4 :

 34 FIGURE 3.4: Model complexity versus bias/variance: blue is the training dataset, orange is the testing dataset. (b) underfitting (blue line, linear fit); (c) reasonable fitting (black curve, third-order polynomial); and (d) overfitting (read curve, fifth-order polynomial).

FIGURE 3

 3 FIGURE 3.5: (a) Confusion matrix of a classification problem with two classes, 0 (negative) and 1 (positive). Each column represents the state in a predicted class, while each row represents the state in an actual class (or observations). Each metric defined here (i.e. TP, FP, FN, TN) can be used to quantify the performance of a classifier. (b) Performance metrics derived from the values in the confusion matrix, including: accuracy, precision, recall, and F1 score.

FIGURE 3 . 6 :

 36 FIGURE 3.6: Illustration of ROC curves: ROC 1 for model 1, ROC 2 for model 2.

  FIGURE 3.7: Folding strategy for the LOO procedure with n f old = n samp folds (corresponding to the n f old rows). In each fold, there are one testing and (n samp -1) training samples.

  FIGURE 3.8: Folding strategy for the LTO procedure with n f old = n samp ×(n samp -1) folds (corresponding to the n f old rows). In each fold, there are one testing, one validation, and (n samp -2) training samples.

  FIGURE 3.9: Illustration of the LTO procedure to estimate a model quality for a fixed complexity level λ with n mod candidate models (horizontal axis). (a) The model errors obtained for each candidate model and each fold of the database B (vertical axes); (b) The obtained RMSE values for the validation and testing datasets. (See detailed description in Sect. 3.7.2.)

  FIGURE 3.10: Schematic illustration of validation and testing RMSE values of predicted yield anomalies for an increasing number of inputs, obtained from the LTO procedure. For a fixed complexity level defined here as the number of inputs, i.e. n input = λ, two RMSE values are obtained: V λ for validation and T λ for testing datasets.

  4.4. This thesis chapter is the subject of two publications in the following international journals: Frontiers in Environmental Science [Dinh et al., 2022b] and Geoscientific Model Development [Dinh and Aires, 2022a].

  FIGURE 4.1: (a) Maps showing averaged coffee area (in 10 3 ha), averaged from 2010 to 2018 over five selected states -with high-intensity Arabica coffee production -in Brazil. (b) Same as (b) but for 20 Robusta coffee-producing districts in the Central Highlands of Vietnam.

Figure 4 .

 4 2 provides two examples showcasing broad phenological stages for coffee in Brazil and the Central Highlands of Vietnam. As a general rule, mature coffee trees (both Arabica and Robusta) will undergo several processes before harvesting: DEC JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC JAN FEB MAR APR MAY JUN JUL AUG ...

.

FIGURE 4

 4 FIGURE 4.2: (a) Climate profile (3-month average of mean temperature and precipitation) in Minas Gerais averaged from 1991 to 2018 and the phenology of Arabica coffee (modified from De Camargo and De Camargo [2001]) corresponding to the tropical climatic conditions of Brazil. (b) Same as (a) but for the tropical monsoon climatic conditions of the Central Highlands in Vietnam, averaged from 2000 to 2018, for Robusta coffee.

FIGURE 4

 4 FIGURE 4.3: (a1) The coffee yield time series and its trend (in 10 3 kg•ha -1 ); and (a2) its corresponding yield anomalies (in 10 3 kg•ha -1 ) of five study states in Brazil, from 1991 to 2018. (b1-2) Same as (a1-2) but for study districts in Dak Lak, from 2000 to 2018.

  FIGURE 4.4: Scheme representing the weather-to-yield impact model used in this section (Sect. 4.2).

FIGURE 4 . 5 :

 45 FIGURE 4.5: Correlation between the yield anomalies and direct weather variables for Minas Gerais (Brazil) and Dak Lak (Vietnam).

FIGURE 4 . 6 :

 46 FIGURE 4.6: Same as Fig. 4.5 but for bioclimatic variables.

FIGURE 4 . 7 :

 47 FIGURE 4.7: The two months with the highest correlation between weather variables (i.e. precipitation and temperature) and coffee yield anomalies, at the district scale, in the Central Highlands of Vietnam: (a) the month with the highest correlation, and (b) the month with the second highest correlation.

  Figure 4.8 shows the final generalization score, i.e. the correlation between

FIGURE 4

 4 FIGURE 4.8: (a) Generalization correlation between observed and estimated yield anomalies; (b) the corresponding explained variance (in %) over selected districts in the Central Highlands (Vietnam). The district-average values are indicated on each panel; AvgDL is the value only over Dak Lak province.

  FIGURE 4.9: (a) The observed (solid lines) and forecast (dashed grey lines) coffee yield anomalies time series and its uncertainties in Cu M'gar (Dak Lak). (b) Same for the coffee yield data (in kg•ha -1 ).

FIGURE 4 .

 4 FIGURE 4.10: The comparison among (1) the true observation, (2) average state from the long-term trend (no weather information), and (3) our weather impact model estimation of the yield anomalies (in 10 3 kg•ha -1 ) and production (in 10 6 kg) over the selected districts in the Central Highlands (Vietnam) in 2011. The absolute values (in 10 6 kg) of the difference between the true production and the estimated production induced from the average state are shown in b3 and from our weather impact model in c3. The district-total values are indicated on panels b2, b3, c2, and c3.

FIGURE 4 .

 4 FIGURE 4.11: Generalization score of different models when applying over different spatial scales (i.e. province and district). The district-average values are indicated on each panel.

FIGURE 4 .

 4 FIGURE 4.12: The training, validation, and/or testing RMSE values of the predicted coffee yield anomalies, using different LIN models (with 18 potential predictors) by increasing the number of inputs, in Cu M'gar (Dak Lak, Vietnam): (a) is induced from LOO procedure, (b) is from LTO procedure.

FIGURE 4 .

 4 FIGURE 4.13: The training, validation, and/or testing RMSE values of the predicted coffee yield anomalies, using LIN5 models by increasing the number of potential predictors, in Cu M'gar (Dak Lak, Vietnam): (a) is induced from LOO procedure, (b) is from LTO procedure.

FIGURE 4 .

 4 FIGURE 4.14: The observed (solid line) and LTO estimated (dashed line) coffee yield anomalies time series in Cu M'gar (Dak Lak, Vietnam).

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.1.2 Data processing . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.1.3 Statistical yield models . . . . . . . . . . . . . . . . . . . . . . 96 5.1.4 Model selection and evaluation . . . . . . . . . . . . . . . . . 96 5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.2.1 Yield model selection -Focus on Bas-Rhin and Landes . . . 96 5.2.2 Reliability model assessment . . . . . . . . . . . . . . . . . . 98 5.2.3 Seasonal yield forecasting . . . . . . . . . . . . . . . . . . . . 101 5.3 Conclusions of the chapter . . . . . . . . . . . . . . . . . . . . . . . . 103

FIGURE 5 . 1 :

 51 FIGURE 5.1: Grain maize database: (a) the average planted area (in 10 3 ha), (b) the average production (in 10 5 tons), and (c) the average yield (in ton•ha -1 ) over 96 French departments (dark lines are regions); (d) same as (c) but presenting only 10 major grain-maize-producing departments. All data are averaged from 2000-2010.

FIGURE 5 . 2 :

 52 FIGURE 5.2: The training, validation, and testing RMSE values of the predicted grain maize yield anomalies, using different LIN models (with 12 potential predictors) by increasing the number of inputs in (a) Bas-Rhin and (b) Landes.

FIGURE 5 . 4 :

 54 FIGURE 5.4: The training, validation, and testing RMSE values of the predicted grain maize yield anomalies, using different models by increasing the number of potential predictors, in Bas-Rhin (France): (a) LIN3 and (b) NN3 (with n neuron = 7) models.

FIGURE 5 . 5 :

 55 FIGURE 5.5: The true testing RMSE maps of predicted grain maize yield anomalies in France for LOO (a1, a2) and LTO (b1, b2) approaches, induced from two LIN models with a different number of inputs: LIN3 (a1, b1) and LIN5 (a2, b2).

FIGURE 5 . 6 :

 56 FIGURE 5.6: Boxplots of residuals (the difference between the observed and estimated yield anomalies) for 10 major grain-maize-producing departments: red horizontal bars are medians, boxes show the 25th-75th percentiles, error bars depict the minimum and maximum values, and red + signs are suspected outliers.

FIGURE 5 . 7 :

 57 FIGURE 5.7: The observed (a obs ) and the estimated yield anomalies time series, for different forecasting months from June to September (e.g. a Jun means June forecasting), for grain maize in Landes (France).

FIGURE 6 . 1 :

 61 FIGURE 6.1: Study area [8 • W -31.5 • E; 33 • N -70 • N] covering Europe and a small part of Africa: (a) the topography elevation (in m), (b) January mean 2 m temperature (in • C), and (c) January total precipitation (in mm) in 2014.

  FIGURE 6.4: Schematic presenting the two main objectives of a calibration method: (1) Do the calibrated data (X C,h ) get closer to the observed data (X O,h )? (2) How the calibrated change (D Cal ) is compared to the raw model change (D M od ). D represents the climate change signals (between historical and future climate), these changes can be measured using several statistical characteristics (e.g. mean, standard deviation, skewness).

Figure 6 .

 6 Figure 6.5 presents an example of July statistics for standard deviation (a) and skewness (b) over Europe and a small part of Africa.Here we can see that the raw climate model data (X M,h and X M,f ) change significantly over time. For example, the July temperature variability increases about 2 • C over the study area (shown in Fig.6.5a). In addition, the skewness varies a lot between the two considered periods. In Fig.6.5b, the skewness of the July temperature series over France is positive for the historical period (panel of X M,h ), i.e. right skewness -most values are on the left side of the mean value. However, this value is mostly negative for the future period (panel of X M,f ). These changes, for instance, in standard deviation and skewness, suggest that a calibration method should better correct

FIGURE 6 .

 6 FIGURE 6.5: July temperature statistics over Europe and a small part of Africa, for standard deviation (a) and skewness (b). For both metrics, we show the observed (Obs), raw model (Mod), and calibrated data (Cal) induced from three quantilebased methods (i.e, FQM, AQM, SDM), for the historical (1981-2014) and future (2066-2099) periods.

Figure 6 .

 6 6 shows the impact of calibration methods on climate change signals D for July temperature, over Europe and a small part of Africa, between the future and historical periods. Most of the calibration methods (except FQM) preserve

FIGURE 6 . 6 :

 66 FIGURE 6.6: Impact of calibration methods on climate change signals D for July temperature between the future (2066-2099) and historical (1981-2014) periods over Europe and a small part of Africa. The first column shows the raw model change (D M od ), other columns present the calibrated changes (D Cal ) induced from several methods. The figure presents the changes in various characteristics of the distribution (mean µ, standard deviation (std) σ, 10 th percentile X10, 90 th percentile X90, skewness γ) in addition to the trend. The spatial correlation between D M od and D Cal are indicated on each panel.

FIGURE 6 . 7 :

 67 FIGURE 6.7: Same as Fig.6.6, but for July precipitation. The coefficient of variation (CV) diagnostic is added, and there is no trend diagnostic for precipitation.

FIGURE 6 .

 6 FIGURE 6.8: Performance of six calibration methods on the climate change signals D between the future (2066-2099) and historical (1981-2014) periods, for temperature (top) and precipitation (bottom). The RMSE between the raw model change (D M od ) and calibrated change (D Cal ) for different moments of the distribution in addition to trends are presented here. This score is computed over all pixels shown in Fig. 6.6 (Europe and a small part of Africa) for monthly calibration. All the months -January (J), February (F), • • • , December (D) -are presented.

FIGURE 6 .

 6 FIGURE 6.10: Schematic of the 3 × 3 moving window to train the calibration model. (1) We use the pixel of interest and its surroundings to build a database that is used to train the local calibration model. (2) The localized transfer model is used to calibrate the central pixel.

FIGURE 6 .

 6 FIGURE 6.11: Comparison of the performance of the pixel-wise and moving window configurations, for January temperature. From the left to the right: the model data (Mod); the EqQM calibrated data (Cal) for the pixel-wise and moving window configurations; the (Cal -Mod) difference and a zoom. Data are averaged for the 2076-2087 period.

FIGURE 7 . 1 :

 71 FIGURE 7.1: Arabica area: (a1) the South America domain with a corresponding study area marked; (b1) the study area (57.5 -37 • W and 12 -28 • S) and topography elevation (in m); (c1) same as (b1), but showing the soil type. Robusta area: (a2, b2, c2) same as (a1, b1, c1), but for the Southeast Asia domain and its corresponding study area (104.22 -109.94 • E and 9.02 -17.82 • N).

FIGURE 7

 7 FIGURE 7.3: Normalised histogram of the suitability prediction, for (a) Arabica and (b) Robusta coffee cases. A threshold of 0.5 is chosen for both cases.

FIGURE 7 . 4 :

 74 FIGURE 7.4: Current and future suitability for Arabica. The future suitabilities are computed from the GERICS_NCC_0 model (i.e. model No. 1 presented in Tab. 7.2) under the RCP2.6 scenario. Dark blue (or 100 %) indicates high suitability, and white (or 0 %) means low suitability. The corresponding suitable areas (in 10 4 ha) are indicated in each panel. The area is computed by summing the entire land area within suitable cells.

FIGURE 7 . 5 :

 75 FIGURE 7.5: Same as Fig. 7.4, but for Robusta. The mean and standard deviation (std) of future suitabilities are computed from eight RCMs (i.e. models No. 2 to 9 presented in Tab. 7.2) under the RCP8.5 scenario.

Figs

  Figs. 7.5b2, c2, and d2. Compared to the current suitability, highly suitable areas will decrease significantly in the future. These areas become smaller and smaller in time. For example, the suitable area is about 81 × 10 4 ha for the current period (Fig.7.5a); however, it will be reduced by half in the 2065s, i.e. ≈ 42 × 10 4 ha (Fig.7.5d1). In addition, the eight simulations give similar predictions by showing small standard deviation values, i.e. 15 % on average, for all three projected periods.

FIGURE 7 . 6 :

 76 FIGURE 7.6: Suitability Arabica change (future suitability -current suitability) for three future periods (2038-2048, 2049-2059, and 2060-2070) in the RCP2.6 and RCP8.5 scenarios. Dark magenta (or -100 %) presents areas with drastic changes from suitable to unsuitable and dark green (or 100 %) indicates positive changes.

FIGURE 7 . 7 :

 77 FIGURE 7.7: Same as Fig. 7.6, but for Robusta.

  FIGURE 7.8: Distribution of suitable Arabica areas by elevation for three future periods (i.e. 2038-2048, 2049-2059, and 2060-2070). Black lines indicate the current (2009-2019) suitable area. Colour lines present the future suitable areas induced from the GERICS_NCC_0 model (i.e. model No. 1 in Tab. 7.2), under RCP2.6 and RCP8.5 scenarios. The areas are calculated by summing the entire land area within suitable cells across and 100 m elevation classes.

FIGURE 7 . 9 :

 79 FIGURE 7.9: Same as Fig. 7.8, but for Robusta. Colour lines and corresponding shaded areas present the mean and standard deviation of future suitable areas induced from six RCMs (i.e. models No. 4 to 9 in Tab. 7.2).

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 8.2 Discussions and perspectives . . . . . . . . . . . . . . . . . . . . . . 154

  Chapter 2 Analysis of available datasets

	No. French name	English name	No. French name	English name
	Blé tendre d'hiver	Winter soft wheat	14 Maïs grain	Grain maize
	Blé tendre de printemps	Spring soft wheat	15 Maïs grain irrigué	Irrigated grain maize
	Total blé tendre	Total soft wheat	16 Maïs grain non irrigué	Non-irrigated grain maize
	Blé dur d'hiver	Winter durum wheat	17 Maïs semence	Maize seed
	Blé dur de printemps	Spring durum wheat	18 Maïs (grain et semence)	Maize (grain and seed)
	Total blé dur	Total durum wheat	19 Sorgho	Sorghum
	Seigle et méteil	Rye and meslin	20 Triticale	Triticale
	Orge et escourgeon d'hiver	Winter barley and bere 21 Autres céréales non mélangées	Other unmixed cereals
	Orge et escourgeon de printemps Spring barley and bere 22 Mélanges de céréales (hors méteil) Mixed cereals (excluding meslin)
	10 Total orge et escourgeon	Total barley and bere	23 Céréales (sauf riz)	Cereals (except rice)
	11 Avoine d'hiver	Winter oats	24 Riz Indica	Indica Rice
	12 Avoine de printemps	Spring oats	25 Riz Japonica et autres riz	Japonica rice and other rice
	13 Total avoine	Total oats	26 Riz	Rice
			27 Toutes céréales	All cereals

TABLE 2 .

 2

2: Available crop data on yields, total production, and harvested areas, at the French regional (department) level. All data are available from 1989 to 2010, except irrigated grain maize (1995-2010) and non-irrigated grain maize

(2007)(2008)(2009)(2010)

.

FIGURE 2.8: Crop yield time series (in t•ha -1 ) of 96 French departments from 1989 to 2010 for winter wheat, winter durum wheat, winter oat and spring oat.

TABLE 2

 2

	Number Crop (French name)	Seasonal type Years with data
	1	Soft wheat (Froment, Blé tendre) Winter	1943-2016
	2		Spring	1943-2016
	3		Total	1900-2016
	4	Durum wheat (Blé dur)	Total	1961-2016
	5	Barley (Orge)	Winter	1943-2016
	6		Spring	1943-2016
	7		Total	1900-2016
	8	Oats (Avoine)	Winter	1943-2016
	9		Spring	1943-2016
	10		Total	1900-2016
	11	Maize (Maïs)	(none)	1900-2016
	12	Potatoes (Pommes de terre)	(none)	1900-2016
	13	Rape (Colza)	Winter	1944-2016
	14		Spring	1943-2016
	15		Total	1900-2016
	16	Sugar beet (Betterave)	(none)	1900-2016
	17	Sunflower (Tournesol)	(none)	1944-2016
	18	Wine (Vigne)	(none)	1900-2016

.3: Dataset description for yields on department level, presented in the study of

[START_REF] Schauberger | Yield trends, variability and stagnation analysis of major crops in France over more than a century[END_REF]

in 2020, preceded only by Asia

[FAO, 2021]

. The major producers of wheat in Europe are Russia, France, Ukraine, and Germany. The Food Security Unit of the Joint Research Center (JRCD5) has just published an exciting dataset: the harmonised sub-national wheat statistics [Food Security Unit of the Joint Research Center (JRC.D.5), 2021]. The sub-national unit refers to different NUTS levels (i.e. NUTS version 2016). This dataset is collected for the EU from National Statistical Institutes and the Eurostat REGIO DB, and it can be downloaded from https://agri4cast.jrc.ec.europa.eu/DataPortal/. Yearly data on area, production, and yield are available from 1975 to 2019 for three categories:

  Thus, we will need to have data that are compatible with what we can use or what is available Map of 13 meteorological stations (red-filled circles) over the Central Highlands of Vietnam, and (b) corresponding locations (longitude, latitude, and height above the mean sea level (i.e. elevation) of these stations. The figure is adopted from the study of Ngo-Thanh et al. [2018].

		Chapter 2 Analysis of available datasets
	(a)	(b)
	Methodology 3.1 Statistical/machine learning models . . . . . . . . . . . . . . . . . . FIGURE 2.14: (a) Chapter 3 3.1.1 Linear regression . . . . . . . . . . . . . . . . . . . . . . . . .
		3.1.2 Random forest . . . . . . . . . . . . . . . . . . . . . . . . . .
		3.1.3 Neural networks . . . . . . . . . . . . . . . . . . . . . . . . .
		3.1.4 Linear mixed-effects . . . . . . . . . . . . . . . . . . . . . . .
		3.2 Model complexity and overfitting . . . . . . . . . . . . . . . . . . . .
		3.2.1 Bias, variance, and model complexity . . . . . . . . . . . . .
		3.2.2 Overfitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
		3.3 Model performance evaluation . . . . . . . . . . . . . . . . . . . . .
		3.4 The training, validation, and testing datasets . . . . . . . . . . . . .
		3.5 Regularization techniques . . . . . . . . . . . . . . . . . . . . . . . .
		3.5.1 Variable selection . . . . . . . . . . . . . . . . . . . . . . . . .
		3.5.2 Loss function modification . . . . . . . . . . . . . . . . . . .
		3.5.3 Early stopping . . . . . . . . . . . . . . . . . . . . . . . . . .
		3.6 Dimensionality reduction . . . . . . . . . . . . . . . . . . . . . . . .
		3.7 Cross-validation (when having limited data) . . . . . . . . . . . . .
		3.7.1 Traditional leave-one-out . . . . . . . . . . . . . . . . . . . .
		3.7.2 Proposed leave-two-out . . . . . . . . . . . . . . . . . . . . .

in a climate change experiment (presented in Sect. 2.2.2). In this case, a model trained by data with high quality (i.e. a precise point measurement) might not be compatible with what is available for climate change applications. Furthermore, the point measurements of a station might not be representative of the pixel value of a model, so the calibration might be difficult.

  n samp = number of samples; %years n pre = number of potential predictors; n mod = number of models; n f old = number of folds of the dataset;

	Score(2,n f old ,n mod ); %representing RMSE or COR; 2 for [Test,Val];
	bm = best model	{1, • • • , n mod };
	%Step 1: Build scores for each fold, each model
	for inp = 1 to n f old
	%Define the folding process
	Test = 1 sample	{1, • • • , n samp };
	Val = 1 sample		{1, • • • , n samp } -Test;
	Learn = {1, • • • , n samp } -Test -Val;
	for imod = 1 to n mod
	%Train models
	model = train(model, Learn);
	Score(1,inp,imod) = RMSE(model, Test);
	Score(2,inp,imod) = RMSE(model, Val);
	end
	end
	%Step 2: Choose best model for all folds; estimate its score
	for isamp = 1 to n

samp M ean V al = mean(Score(2,n f old {isamp},:)); %(1,1,n mod ) ibm(isamp) = argmin i (M ean V al ); Score T est (isamp) = mean(Score(1,n f old {isamp},ibm(isamp))); Test score Score V al (isamp) = mean(Score(2,n f old {isamp},ibm(isamp))); Val score end F inalScore T est = mean(Score T est ) F inalScore V al = mean(Score V al )

  Discussions and conclusions of the chapter . . . . . . . . . . . . . . 90

4.1 Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.1.1 Study area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.1.2 Overview of coffee phenology . . . . . . . . . . . . . . . . . 65 4.1.3 Data processing . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2 Weather-to-yield impact model using LOO and regularization techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.2.1 Weather-to-yield impact model . . . . . . . . . . . . . . . . . 72 4.2.2 Correlation analysis . . . . . . . . . . . . . . . . . . . . . . . 73 4.2.3 Results of the impact model . . . . . . . . . . . . . . . . . . . 77 4.2.4 Yield estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.2.5 Sensitivity to the model spatial scale . . . . . . . . . . . . . . 82 4.2.6 Discussions of the Robusta analysis . . . . . . . . . . . . . . 84 4.3 Weather-to-yield impact model using LTO . . . . . . . . . . . . . . . 86 4.3.1 Statistical yield models . . . . . . . . . . . . . . . . . . . . . . 86 4.3.2 Yield model selection . . . . . . . . . . . . . . . . . . . . . . . 87 4.3.3 Yield anomaly estimation . . . . . . . . . . . . . . . . . . . . 89 4.4

Bud development Flowering Fruit development Maturation ... Bud formation Dormacy Pinhead Rest ... Bud growth Grain expansion Grain formation Vegetative growth Maturation Harvest Productive stage Vegetative growth for the next harvest ... Timeline corresponding to the tropical climatic conditions of Brazil ... H -21 H -20 H -19 H -18 H -17 H -16 H -15 H -14 H -13 H -12 H -11 H -10 H -9 H -8 H -7 H -6 H -5 H -4 H -3 H -2 H -1 H 0 (a) T (in °C) P (in mm) ... MAR APR MAY JUN JUL AUG SEP OCT NOV DEC JAN FEB MAR APR MAY JUN JUL AUG SEP OCT NOV DEC JAN FEB …

	..	Y(t-1)	Y(t)	...

... Y(t-1) Y(t) Y(t+1) Irrigation ... ..

. Bud development Flowering Fruit development Maturation Bud formation Dormancy Pinhead Rest Bud growth Dormancy Grain expansion Grain formation Vegetative growth Maturation Harvesting Productive stage Vegetative growth for the next harvest … Timeline corresponding to the tropical monsoon climatic conditions of the Central Highlands (Vietnam)

	(b)
	…	H -21	H -20	H -19	H -18	H -17	H -16	H -15	H -14	H -13	H -12	H -11	H -10	H -9	H -8	H -7	H -6	H -5	H -4	H -3	H -2	H -1	H 0

TABLE 4

 4

.1: Bioclimatic variables which corresponds to the period from the bud development process to the harvest season's peak (as shown in Fig.

4

.2). Here, n = 19 results in 38 potential inputs (i.e. 19 realisations (chronological data) for P and T series), for modelling coffee yield.

TABLE 4

 4

.2: Correlation (COR) between observed and estimated yield anomalies over the Dak Lak districts.

TABLE 5 .

 5

1: The correlation between the observed and estimated yield anomalies for different forecasting months (from June to September), over 10 major grainmaize-producing departments.

  [START_REF] Cannon | Bias correction of gcm precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes[END_REF] DetQM better preserves the model's long-term trend in the mean.EqQM tries to estimate the correction on the historical period but only on the subspace of the variability that is predictable by a climate model.[START_REF] Li | Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change ar4 models using equidistant quantile matching[END_REF][START_REF] Wang | Equiratio cumulative distribution function matching as an improvement to the equidistant approach in bias correction of precipitation[END_REF][START_REF] Pierce | Improved bias correction techniques for hydrological simulations of climate change[END_REF] Using the QM that links the future climate to the historical, EqQM can isolate only the longterm tendency from the climate model, but not the other variabilities. In contrast, e.g. for the direct QM, the calibration may take into account many variabilities from observations (both temporal and spatial) that we do not expect to have in the model.

	Equidistant/Equiratio
	Quantile Mapping
	(EqQM)
	Quantile Delta Mapping (QDM)	QDM shares the same ideas/results as EqQM.	Cannon et al., 2015
		SDM is conceptually similar to QDM (or EqQM) but more explicitly considers the likelihood
	Scaled Distribution Mapping (SDM)	of events (i.e. recurrent interval). Also, SDM includes the preprocessing step that deals with the rain-day frequencies and trends in temperature values. SDM focuses on preserving the raw climate model change in terms of magnitude, rain-day	Switanek et al., 2017; Maraun et al., 2021
		frequencies, and the likelihood of events.

TABLE 6 .

 6 1: Summary of the calibration methods.

  AI and, more specifically, machine learning techniques are vastly spreading in the Geoscience discipline. These types of data-driven techniques could potentially be used to calibrate climate simulations. In fact, AI and traditional techniques are regression models based on a learning database, equivalent to a calibration transfer function set up on the historical record. For instance, neural networks (NN) models have the advantage of being non-linear models, which could be helpful to represent quantile-based transfer functions, potentially non-parametric, with no required hypothesis in the distributions. This NN representation richness and flexibility can be an advantage, but it is dangerous too because what has been done in the calibration is less controlled and understood. The non-linearity of AI models also means that the NN transformation would be state-dependent, which could be an exciting feature of the calibration. Furthermore, NN deals extremely well with multivariate relationships, but again this can be risky as inferring multivariate relationships from data can be driven by most representative cases, at the detriment of the rarest cases. This could conduct the calibration to smooth out the climate simulations and simplify multivariate relationships among climate variables. Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 7.1.1 Study area and coffee data . . . . . . . . . . . . . . . . . . . . 138 7.1.2 Potential predictors . . . . . . . . . . . . . . . . . . . . . . . . 140 7.1.3 Climate data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 7.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.2.1 Learning database identification . . . . . . . . . . . . . . . . 141 7.2.2 Suitability model . . . . . . . . . . . . . . . . . . . . . . . . . 142 7.2.3 Calibration of climate simulations . . . . . . . . . . . . . . . 145 7.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 7.3.1 Model validation and variable contribution . . . . . . . . . . 145 7.3.2 Climate suitability . . . . . . . . . . . . . . . . . . . . . . . . 146 7.4 Discussions and conclusions of the chapter . . . . . . . . . . . . . . 151

	Chapter 7
	Climate change impacts on coffee
	production in Brazil and Vietnam
	7.1

  Chapter 7 Climate change impacts on coffee production in Brazil and Vietnam

	Study case No.	Abbreviation of RCM experiment	Driving model (GCM) Realisation RCM	RCP2.6 RCP8.5
	Arabica	1	GERICS_NCC_0	NCC-NorESM1-M (NCC, Norway)	r1i1p1	REMO2015 (GERICS, Germany)	x	x
		2	SMHI_CNRM	CNRM-CM5 (CNRM, France)	r1i1p1	RCA4 (SMHI, Sweden)		x
		3	SMHI_Had	HadGEM2-ES	r1i1p1	RCA4 (SMHI, Sweden)		x
		4	ICTP_Had	(Hadley Centre, UK)	r1i1p1	RegCM4 (ICTP, Italy)	x	x
	Robusta	5	GERICS_Had		r1i1p1	REMO2015 (GERICS, Germany)	x	x
		6 7	ICTP_NCC GERICS_NCC	NCC-NorESM1-M (NCC, Norway)	r1i1p1 r1i1p1	RegCM4 (ICTP, Italy) REMO2015 (GERICS, Germany)	x x	x x
		8 9	ICTP_MPI GERICS_MPI	MPI-ESM-MR (MPI-M, Germany)	r1i1p1 r1i1p1	RegCM4 (ICTP, Italy) REMO2015 (GERICS, Germany)	x x	x x
	TABLE 7.2: List of regional climate models used in this chapter, for Arabica and
				Robusta cases.

Also known as mechanistic (dynamic), biophysical, or crop simulation models

Here, we refer to the metropolitan area of France -which is geographically in Europe -and Corsica.

at the time we gathered the data for this thesis

at the time of writing this thesis

We often consider the yield anomaly values.

Also called the predictors in the following.

One sample corresponds to one year of yield data.

Also known as artificial neural networks, or neural nets

Also known as mixed models or mixed error-component models.

In general, many statistical tools do folding by default; however, the validation score is rarely used in a rigorous way in the model selection tasks.

Also known as the bias correction[Ho et al.,

2012;[START_REF] Hawkins | Calibration and bias correction of climate projections for crop modelling: An idealised case study over Europe[END_REF] or the nudging approach[START_REF] Hawkins | Calibration and bias correction of climate projections for crop modelling: An idealised case study over Europe[END_REF]].2 Also known as the change factor[START_REF] Diaz-Nieto | A comparison of statistical downscaling and climate change factor methods: impacts on low flows in the river thames, united kingdom[END_REF] Chen et al., 2011a;[START_REF] Ho | Calibration strategies a source of additional uncertainty in climate change projections[END_REF][START_REF] Hawkins | Calibration and bias correction of climate projections for crop modelling: An idealised case study over Europe[END_REF] or the delta change[START_REF] Hay | A comparison of delta change and downscaled gcm scenarios for three mountain basins in the united states[END_REF][START_REF] Maraun | Bias correcting climate change simulations -a critical review[END_REF].

Also known as "probability mapping"[START_REF] Ines | Bias correction of daily gcm rainfall for crop simulation studies[END_REF][START_REF] Block | A streamflow forecasting framework using multiple climate and hydrological models1[END_REF], "quantilequantile mapping/quantile-based mapping"[START_REF] Wood | Hydrologic implications of dynamical and statistical approaches to downscaling climate model outputs[END_REF][START_REF] Déqué | Frequency of precipitation and temperature extremes over France in an anthropogenic scenario: Model results and statistical correction according to observed values[END_REF][START_REF] Boé | Statistical and dynamical downscaling of the seine basin climate for hydro-meteorological studies[END_REF], "statistical bias correction"[START_REF] Piani | Statistical bias correction for daily precipitation in regional climate models over Europe[END_REF][START_REF] Haerter | Climate model bias correction and the role of timescales[END_REF], or "histogram equalization"[START_REF] Rojas | Improving pan-european hydrological simulation of extreme events through statistical bias correction of rcm-driven climate simulations[END_REF].

Chapter 6

Calibration of the climate simulations

applications; it can be the mean, variance, frequency, quantiles, or other complex properties. Moreover, several models can be used to calibrate these properties, from linear scaling to non-linear transformations. This chosen model will be referred to as the transfer model or "transfer function" in the following.

(a) Linear scaling (LS)

The LS method [START_REF] Lenderink | Estimates of future discharges of the river rhine using two scenario methodologies: direct versus delta approach[END_REF]Graham et al., 2007a;[START_REF] Berg | Bias correction of high resolution regional climate model data[END_REF]] is one of the simplest calibration methods which aims to preserve the mean value of the considered variable: the mean of calibrated model data will be equal to that of the observations. A linear transformation model is used. An additive (or multiplicative) term corresponding to the difference of historical mean observations and historical model simulations is applied to the model temperatures (or precipitations).

The calibrated temperatures can be described as:

Chapter 7 Climate change impacts on coffee production in Brazil and Vietnam

Potential predictors

We chose a set of 12 variables (Tab. 7.1) as the potential model predictors. These potential predictors are derived from precipitation, temperature, and evaporation, characterising growing conditions, e.g. between flowering (i.e. from January to March), growing season (i.e. from March to September), or harvest period (i.e. from October to December). They have been selected according to expert knowledge and recommendations [Bunn et al., 2015b;[START_REF] Lambot | Chapter 2 -cultivating coffee quality-terroir and agro-ecosystem[END_REF][START_REF] Kath | Not so robust: Robusta coffee production is highly sensitive to temperature[END_REF][START_REF] Kath | Temperature and rainfall impacts on robusta coffee bean characteristics[END_REF]Dinh et al., 2022b].

Climate data

Monthly data on total precipitation, mean temperature, maximum and minimum daily temperature, and evaporation are considered for both current and future climate. For the current climate , we used the ERA5-Land reanalysis dataset (Sect. 2.2.1). We then integrated these data from its original 0.1 • × 0.1 • resolution into 0.22 • × 0.22 • to be consistent with the future climate data, which will be presented in the following.

Regional climate simulations and projections are obtained from several model simulations of the CORDEX-SAM22 and CORDEX-SEA22, as presented in Sect. 2.2.2. The list of regional climate models (RCMs) carried out in this chapter is presented in Tab. 7.2. More models should have been included, especially for the Arabica study; however, we only test here nine models due to the time limitations of the thesis. We considered the historical period from 1981 to 2005, used for the calibration in Sect. 7.2.3. We will investigate the future period from 2031-2077 for the RCP scenarios of the CMIP5: the low (RCP2.6) and high (RCP8.5) greenhouse gas concentration scenarios [START_REF] Van Vuuren | The representative concentration pathways: an overview[END_REF]. accounts for four inputs. With a higher number of inputs, the model is overfitted by poor generalization ability [START_REF] Dinh | Climate change impacts on robusta coffee production in vietnam[END_REF].

Model training and evaluation

To assess the model's generalization, we divided our database into three sets: training (60 %), validation (20 %), and testing (20 %). The model is then evaluated using performance metrics derived from the confusion matrix (see more in Sect. 3.3), which are commonly used for the evaluation of classification models. We considered here two common metrics:

• precision: the measure of correctly identified coffee samples over the number of all correctly identified samples;

• recall: the measure of correctly identified coffee samples over the count of actual coffee samples.