Viticulture extension in response to global climate change drivers – lessons from the past and future projections

The potential areal extent of agricultural crops is sensitive to climate change and its underlying drivers. To distinguish between the drivers of past variations in the Mediterranean viticulture extension since Early Antiquity and improve projections for the future, we propose an original attribution method based on an emulation of offline coupled climate and ecosystem models. The emulator connects the potential productivity of grapevines to global direct and indirect climate drivers, notably orbital parameters, solar and volcanic activities, demography, and greenhouse gas concentrations. This approach is particularly useful to place the evolution of future agrosystems in the context of their past variations. We found that variations in potential area for viticulture during the last 3 millennia in the Mediterranean Basin were mainly due to volcanic activity, while the effects of solar activity and orbital changes were negligible. In the future, as expected, the dominating factor is the increase in greenhouse gases, causing significantly drier conditions and thus major difficulties for viticulture in Spain and North Africa. These constraints will concern significant areas of the southern Mediterranean Basin when global warming exceeds +2 ^∘C above preindustrial conditions. Our experiments showed that even intense volcanic activity comparable to that of the Samalas – sometimes considered to be the starting point of the Little Ice Age in the mid-13th century – would not decrease aridity and so not slow down this decline in viticulture extension in the southern margin of the Mediterranean area. This result does not confirm the idea of geoengineering that solar radiation modification (SRM) is an efficient option to limit future global warming.

1 Introduction

Our scientific question is related to the attribution of viticulture, which has had an economic role in the Mediterranean Basin since antiquity, extension changes to any natural or anthropogenic drivers. The cultivation and domestication of the grapevine began between the 7th and 4th millennia before the common era (BCE) between the eastern Mediterranean and Caspian areas and spread to Egypt, the Middle East, and the entire Mediterranean (Terral et al., 2010; Bouby et al., 2021). Introduced in the Gaul region (i.e., roughly France and surrounding regions) by Greek colonists ca. 600 BCE, around the time they settled Marseilles, viticulture was initially limited to Mediterranean Gaul (Bouby et al., 2014). Vineyards expanded into the northern part of Gaul in the 1st century CE, where wine production developed quite considerably in the following centuries up to the Paris region, the Rhine and Moselle valleys (Brun, 2010), and even in southern England (Brown et al., 2001). One hypothesis behind this expansion is the climate warming during the Roman Climatic Optimum (RCO) (Mccormick et al., 2012). These climate variations are driven by global forcing variables such as solar or volcanic activity (Wanner et al., 2008; Brayshaw et al., 2010; Fuks et al., 2017). After the RCO, the temperature decreased significantly and Gaul entered the so-called Late Antique Little Ice Age, or LALIA (536–660 CE) (Büntgen et al., 2016). This change may have been triggered by several large volcanic eruptions at 536, 540, and 547 CE (Sigl et al., 2015). This assumption remains difficult to prove because of the limited historical and archeological sources. In any case, the 500 to 900 CE period remained relatively cold with oscillating precipitation changes in the region (Reale and Dirmeyer, 2000). In Europe, the following Medieval Climate Anomaly (MCA, approx. 900–1200 CE) (Luterbacher et al., 2016) was likely of intensity comparable to the RCO. The Little Ice Age (LIA, approx. 1250–1850 CE) was a period of alpine glacier advance (Holzhauser et al., 2005), again marked by several large volcanic eruptions, in particular that of the Samalas (Indonesia) in 1257 (Lavigne et al., 2013).

On timescales of centuries and millennia, quality and yield of agricultural crops are strongly affected by climate fluctuations. The nature of the change depends on external forcings and internal feedbacks of the climate system, which produce different spatial and seasonal patterns of the main variables in the atmospheric environment. The main objective of this paper is to develop an innovative solution to statistically model the impact of changing climate forcing on vegetation over several millennia using the grapevine (a major crop of the Mediterranean and European region) as an example. We mimic this impact model on the basis of a large ensemble of existing model simulations using statistical relationships much faster to compute (Kennedy and O'Hagan, 2000). Based on a large range of climate states from high-resolution simulations with coupled Earth system models for the last glacial period to future global warming scenarios, this approach, called a “emulator”, provides robust results and can be applied to a large range of ecosystem processes under different conditions.

The global drivers of the studied changes are anthropogenic, including greenhouse gas (GHG) emissions, land use and cover changes, population density, and economic production, and natural, including the Earth's orbit as well as volcanic and solar activity. They are the boundary conditions of Earth system models (ESM) (Kay et al., 2014). The ecosystem processes are assessed using impact models (IMs) coupled to these ESMs (Franklin et al., 2016; Warszawski et al., 2014; Frieler et al., 2017). In most cases, coupling is offline because (i) the spatial scales of ecosystem processes are much finer than are those of ESMs, thus making a scale transfer necessary (Su et al., 2016); (ii) a climate simulation of a given ESM is interpreted as one realization out of a set of possibilities determined by the boundary conditions and the characteristics of the ESM, thus making it necessary to work on ensembles of models to be representative of the climate system (Kay et al., 2014); and (iii) each ESM has intrinsic biases that must be corrected before it can be used to drive the ecosystem model (Williamson et al., 2015).

Even if human innovation and colonization were also responsible for the expansion of viticulture, the climate must be suitable for grape cultivation and thus remains a control variable. As soon as the climate changes to worse conditions for wine growth, it becomes a driver of a decline in viticulture. Such fluctuations are particularly noticeable near the northern range limit of wine growth during the end of the Roman and Medieval periods with a regression of the grape cultivation.

Although we focus our viticultural analysis on the Gaul region (France and surrounding countries), we need to enlarge the area to the entire Mediterranean and European surrounding region to robustly capture the relationships between global drivers and viticulture extension. For the same reason, we use a large diversity of time slices of the past (Last Glacial Maximum, mid-Holocene, last millennium) and of the future up to 2100 according to several scenarios. The widely diverse situations used for calibration made it possible to produce a robust emulator that was effective for extrapolating a wide range of past and future climate states. This method is appropriate to establish a bridge between past variations of a given agrosystem and its future evolution. Once calibrated, the emulator is of great help to efficiently support decisions in the domain of impacts of climatic change.

2 Material and methods

As ESMs are an imperfect representation of reality, our approach has the same limitations even if the model parameters are carefully chosen (Crucifix, 2012). Crucifix (2012) developed the idea that a Bayesian framework is adequate to stimulate the modeling of the distance between the ESM and reality. So ESMs produce a certain vision of the world that needs to converge to the reality represented by the paleodata. This procedure is called data assimilation, which provides a way to reduce the uncertainties.

The concept of our method and the data used are presented in Fig. 1. Global forcings are the inputs of low-resolution ESMs. Step 1 involves adapting output fields to a common high-resolution grid using statistical downscaling. This step is also useful to correct the systematic biases of the ESM simulations. Step 2 involves applying the vegetation model BIOME4 to the high-resolution fields to calculate ecosystem variables. Steps 1 and 2 are repeated for each ESM simulation available. These simulations based on different forcings and different models provide a large and diversified calibration dataset. They are a guarantee that the emulator is in some way more robust than the ESM simulations taken separately. Calibration of the emulator (step 3) is done by spatial regression of the ecosystem variables on the forcing variables. Step 4 involves the application of the emulator to the forcings of new time slices or future scenarios. During this step, the annual temperature and precipitation results for the past time slice were compared to paleodata to identify the amplitude of the volcanic and solar activity effects (data assimilation). Step 5 is the independent validation of the emulator using tree-ring data.

Figure 1 Diagram of the different steps in the emulator approach; the ecosystem variables, related to bioclimate and net primary productivity, are defined in Table 1. Low-density hashes in the trapezes indicate low resolution of the corresponding data, and high-density hashes correspond to higher resolution.

Table 1 The biome types simulated by BIOME4 and output variables used for the paper; note that there is no correspondence between the two columns.

Figure 2 Natural (volcanic activity and solar activity) and anthropogenic forcing (CO₂) for the last millennium: (a) the effective aerosol radius, which is directly related to the aerosol optical depth, is calculated by Crowley and Unterman (2013); (b) the proxy indicating the total solar irradiance (TSI), i.e., the solar modulation function in MeV based on ¹⁴C (Muscheler et al., 2007); and the (c) atmospheric CO₂ concentration (in ppm). The red horizontal lines represent the mean values of the five selected periods.

2.1 The forcing variables

Climate forcings or drivers are perturbations imposed on the Earth's energy balance. They are mainly the following.

• i.

  Orbital parameters: eccentricity (ecc) and obliquity of the ecliptic (obl) as well as solar longitude at the perihelion (ω), or more precisely the climate precession pr = ecc × sin (ω). Nevertheless, to avoid redundancy, we use as forcing variables the three basic parameters (ecc, obl, ω). They drive the orbit of the Earth and influence its climate by modulating the solar energy intercepted by the planet and its seasonal distribution at timescales above thousands of years (Berger and Loutre, 1991). They had a major impact on the Earth's climate from the last glacial period to the Holocene. At the century timescale, the orbital parameters do not operate significantly, so the values of the last millennium can be linearly interpolated between the values at 1000 yr BP and the present values, and the values of the 21st century can be set to the present values.

• ii.

  Greenhouse gas (GHG) concentration: carbon dioxide (CO₂), methane (CH₄), and nitrous oxide (N₂O). The effect of these GHGs has been significant from the glacial to interglacial periods and has a major effect for the current century (see, for example, Fig. 2 for CO₂).

• iii.

  The world population values taken from the Hyde 3.1 database for the past periods from Klein Goldewijk et al. (2011) and for the future periods from van Vuuren et al. (2011). The population varied from less than 10⁶ at the Last Glacial Maximum (LGM) to 7×10⁹ in 2010 and is projected to be between 9×10⁹ and 12.3×10⁹ in 2100 according to the scenario.

• iv.

  Volcano activity (V) represented by the effective aerosol radius deduced from the aerosol optical depth from ice core sulfate records from both polar regions for the last millennium (Crowley and Unterman, 2013) (Fig. 2). Its value is 0.2 when there is no eruption. The maximum value (0.8) was found in 1258, the year after the Samalas eruption (Lavigne et al., 2013). The 21 ka, 6 ka, and future values were set to the preindustrial values.

• v.

  Solar activity (S) is inferred from ¹⁴C and ¹⁰Be records, the proxy for the total solar irradiance (TSI) (Muscheler et al., 2007) for the last millennium, and varies from 0 to 1200 MeV (Fig. 2). The 21 ka, 6 ka, and future values were set to the preindustrial values.

2.2 The past time climate data (used for calibration)

The past time periods retained for calibration are (1) the Last Glacial Maximum (21 ka) characterized by very low greenhouse gases (GHGs), low ecc, and obl and a pr close to zero determined by a low ecc; (2) the mid-Holocene (6 ka) characterized by GHG close to that of the preindustrial period, low ecc and pr, and a large obl increasing the seasonality of the climate; (3) five periods of the last millennium (Fig. 2) representing various combinations of S and V forcings, including V− S+ (1200–1220), V+ −S (1255–1320), V− S+ (1380–1420), V− $S++$ (1720–1780), and V+ S+ (1810–1840); and (4) the preindustrial (PI) reference period with intermediate GHG concentrations, low ecc, medium obl, medium pr, and large V and S.

The simulations were archived by the Paleoclimate Modelling Intercomparison Project (PMIP), which produced snapshot simulations of the mid-Holocene (6 ka BP) and Last Glacial Maximum (21 ka BP) time slices to study the effects of insolation, greenhouse gases, and massive ice sheets on the climate (https://pmip3.lsce.ipsl.fr/overview/, last access: 15 June 2023) (Braconnot et al., 2012). The models used for intercomparisons were coupled atmosphere–ocean–vegetation models with various levels of complexity and resolution (Table 2).

Table 2 CMIP5 and PMIP3 simulations used in the paper. The first column is the code of model, and the second indicates the institute which built the model; the other columns provide the availability of the scenarios (RCP) or time slices. Hist.: historical period (last millennium), MidHol: mid-Holocene (6 kyr BP), LGM: Last Glacial Maximum (21 kyr BP).

2.3 The future climate data (used for calibration)

The future periods are determined by increasingly high GHG concentrations (depending on the scenario, see below), and the other forcing variables were set to the PI value. The simulations were archived by the Coupled Model Intercomparison Project (CMIP). We used simulations of CMIP5 (Taylor et al., 2012). These include both (i) century-scale integrations, which usually start from a preindustrial control, and climate predictions until the end of the 21st century, as well as (ii) near-term integrations for the next 10 to 30 years, which are initialized using observed ocean and sea-ice conditions. The CMIP5 climate change projections are driven by emission scenarios divided into four classes referred to as Representative Concentration Pathways (RCPs) (Moss et al., 2010). The high-emissions scenario, named RCP8.5 or “business as usual”, refers to a high radiative forcing of emissions at the end of the 21st century (8.5 W m⁻²). The low-emission scenario, which roughly represents that of the Paris Agreement, reaches its maximum value near the middle of the century before decreasing to a level of 2.6 W m⁻² at the end of the century. We only used the intermediate scenario RCP4.5, which refers to a radiative forcing maintained at values of 4.5 W m⁻² at the end of the 21st century. We used 10 time slices interpolated at a 10-year steps (2010, 2020, …, 2100). The forcing variables for the future scenarios were obtained from the website http://tntcat.iiasa.ac.at/RcpDb/ (last access: 15 June 2023).

2.4 The application data

We used our emulator to analyze the response of key ecosystem variables to global forcing. We focused on past periods that were marked by important climate and societal changes (Table 3). The present time slice was defined by the mean values for the 1961–1990 period. For the future, instead of using time slices, we defined the scenarios according to the global temperature signal simulated in the different models using the relationship between global warming and CO₂ concentration (Guiot and Cramer, 2016).

• i.

  The +1.5 ^∘C global warming recommended by the Paris Agreement is reached with a CO₂ concentration of 440 ppm; Fig. 3 shows that the average model simulation reaches this value under scenario RCP2.6 in approximately 2040.

• ii.

  The +2 ^∘C global warming is reached with a CO₂ concentration of 480 ppm and is reached under scenario RCP4.5 in approximately 2050 (Fig. 3).

• iii.

  The +3 ^∘C global warming is reached with a CO₂ concentration of 600 ppm and is reached under scenario RCP8.5 in approximately 2060 (Fig. 3).

• iv.

  The +5 ^∘C global warming is reached with a CO₂ concentration of 900 ppm and is reached under scenario RCP8.5 in about 2090 (Fig. 3).

We obtained the corresponding CH₄ and N₂O concentrations and population sizes from the boundary condition database of CMIP5. The other future forcings are set to the present values.

Table 3 Typical past periods used to calibrate our emulator. The paleoclimatic and societal information is found in the literature given in the references from the late Holocene to the present.

Figure 3 Global annual temperature anomalies (T) simulated by the various models of PMIP3 and CMIP5, according to the reference period average (1961–1990). The left graphic represents the past simulations, and the right graphic represents the present and future projections for three scenarios (RCP2.6 in green, RCP4.5 in orange, RCP8.5 in red). Each dot represents one model simulation. The vertical scales are shifted for better readability. The reference is given by the New et al. (2002) dataset calculated on the 1961–1990 period. The anomalies versus the preindustrial period are obtained by adding approximately 0.5 ^∘C to the 1961–1990 anomalies.

Following Giorgi and Lionello (2008), the study area was divided into nine grid boxes (Fig. 4). We added a 10th box corresponding to the Gallia Narbonensis Roman province (south of France), the key area for the introduction of viticulture in Gaul.

Figure 4 Map of the 10 boxes defined in the Mediterranean region.

We also considered two other scenarios based on boundary conditions of +5 ^∘C except for volcanic and solar activity. The idea is to assess whether volcanic and solar conditions typical of the Little Ice Age can moderate the effect of the strong increase in GHG concentrations. The first additional scenario (labeled +5 CV+) was assigned very substantial volcanic activity (0.8), the highest value in the last millennium, and low solar activity (i.e., 100). In contrast, the second additional scenario (labeled +5 CV−) was assigned low volcanic activity (0.1) and high solar activity (700), corresponding to those of the MCA. The forcing values are listed in Table 2.

2.5 Downscaling (step 1)

ESM resolution is too coarse (>100 km depending on the model) to assess the impact and risk of climate changes on ecosystems or human systems. Different approaches can be used to derive higher-resolution information. We use a statistical downscaling (SD) based on statistical relationships that link large-scale atmospheric variables with local and regional climate variables and that are applied to coarser-resolution model simulations (Su et al., 2016; Grouillet et al., 2016; Levavasseur et al., 2011), thereby providing higher-resolution estimates of climate variables. SD methods are among the less computationally expensive downscaling techniques.

The SD that we use assumes the hypothesis that the fields of climate anomalies do not depend on the grid resolution. So, they can be interpolated at a finer resolution than the initial resolution. We use a common 0.5^∘ high-resolution (HR) grid for each climate simulation. It is provided by the high-resolution dataset of surface climate over global land areas (New et al., 2002) (which has a 10 min resolution and was aggregated at a 30 min resolution). This climatology is based on the 1961–1990 period. Each low-resolution (LR) anomaly field provided by the ESM is downscaled to the HR grid by the bilinear interpolator of platform R (function interp. surface of package fields). This HR anomaly field is added to the 1961–1990 climatological HR field of New et al. (2002). As the simulations are corrected by the differences between control simulations and observations, the main systematic biases of the model are removed. For the PMIP3 simulations (LGM, MidHol, and historical), the control is the preindustrial simulation (PI), so we use the 1900–1931 climatology very close to the PI (1861–1890) used by the recent IPCC reports (Allen et al., 2018). Figure 3 presents the global mean temperature for each period and scenario considered.

2.6 Vegetation modeling (BIOME4, step 2)

The ecosystem model is coupled offline to downscaled climate outputs and provides indications of land vegetation structure and productivity. Here, we used a process-based equilibrium vegetation model, BIOME4 (Kaplan et al., 2002), which has been successfully applied to similar questions before (Guiot and Cramer, 2016). While a wide range of climate simulations are necessary to ensure robust calibration, the use of a single ecosystem model is required to ensure that the same ecosystem variables are calculated for all simulations.

BIOME4 simulates some of the most important processes related to vegetation at the scale of the ecosystem, specifically photosynthesis, respiration, and evapotranspiration. It predicts structure and productivity of broad-scale land ecosystems from monthly temperature and rainfall values as well as annual CO₂ concentration. It uses a photosynthesis scheme that simulates acclimation of plants to changed atmospheric CO₂ by optimization of nitrogen allocation to foliage and by accounting for the effects of CO₂ on net assimilation, stomatal conductance, leaf area index (LAI), and ecosystem water balance. BIOME4 is based on a sufficiently simple description of ecophysiological processes to allow broad-scale application. It represents substantial advantages over niche models because it has not been tuned to reproduce present-day potential vegetation, but rather to correctly simulate the main processes underlying the potential vegetation distribution, which are assumed to have been similar throughout the Holocene. BIOME4 does not account for human land use.

The inputs of the model are monthly average values of temperature, precipitation and sunshine percentage, atmospheric CO₂ concentration, and soil texture. The temperature and precipitation variables are directly obtained from the simulations downscaled at the 0.5^∘ grid. The sunshine percentages are obtained by linear regression on temperature and precipitation (Guiot et al., 2000). The absolute minimum temperature is derived from the mean temperature of the coldest month according to the quadratic equation of Prentice et al. (1992). The soil data are provided by the FAO (Zobler, 1986). The CO₂ atmospheric concentration values are those used by CMIP5 and PMIP3 as boundary conditions for their simulations.

The outputs include net primary production (NPP) of each of the potentially occurring 13 plant functional types (PFTs), the total NPP, and the corresponding “biome type” from a set of 28 broad categories (Table 1). The model also provides several bioclimatic variables (Table 1).

2.7 The calibration of the emulator (step 3)

The numerous model outputs available in the CMIP and PMIP databases (Table 2) are used to calibrate robust statistical approximations of coupled ESMs and BIOME4, called the emulator (Kennedy and O'Hagan, 2000). Various emulators are used in climate science (Tran et al., 2016; Zhu et al., 2015; Rougier and Goldstein, 2014; Castruccio et al., 2014), including in paleoclimatology (Strassmann and Joos, 2018; Joos et al., 1996; Bounceur et al., 2015). Our approach is the first ESM-independent emulator because it is calibrated using a large set of model simulations under very different scenarios. The calibration is based on a geographically weighted regression (Brunsdon et al., 1998), whereby the ecosystem variables are expressed as functions of the global forcing variables.

For any choice of input q vector x (forcings), the climate simulator is $y\left(\mathit{x}s\right)=\left(y_{\mathrm{1}}\right(\mathit{x},s)$, …, y_m(x,s)), where m is the number of (bio)climatic variables whose response is to be analyzed and s is the location where the climate variables must be estimated. There are deterministic and possibly nonlinear functions. This model is a statistical representation of the ESM + BIOME4 model, with the very interesting properties to run very fast and provide an uncertainty description for the whole plausible input space (i.e., conform to physics, internally consistent and reasonable; Amara, 1991):

$$\begin{array}{}\text{(1)}& y_j(\mathit{x},s)=\mathbf{X}{\mathit{\beta }}_j\left(s\right)+e_j\left(s\right)j=\mathrm{1},\mathrm{\cdots },m;s=\mathrm{1},\mathrm{\dots },n,\end{array}$$

where s indicates the grid point out of a total of n points, y_j(x,s) is the anomaly of the climatic variable j (out of m) at location s, i.e., the climate variable at time t minus the climate variable at time 0 (present), X is a matrix with v specified columns of the input global forcings (they do not depend on the location s), and e(s) is a stationary Gaussian variable $N(\mathrm{0},{\mathit{\sigma }}_j^{\mathrm{2}})$. This can be considered a least-square problem wherein the coefficients β_j(s) are estimated to minimize the sum of the squared errors. These coefficients can be written as

$$\begin{array}{}\text{(2)}& \widehat{{\mathit{\beta }}_j}\left(s\right)={\left({\mathbf{X}}^T\mathbf{W}\left(s\right)\mathbf{X}\right)}^{-\mathrm{1}}{\mathbf{X}}^{T}W\left(s\right)y_j,\end{array}$$

with W(s) as a matrix of weights specific to location s such that observations nearer to s are given greater weight than observations further away, and T is the transpose symbol. The estimation of β_j(s) in location s is provided by weighting all the observations according to their distance from s, d_is. The method is called geographically weighted regression (GWR) (Brunsdon et al., 1998). We use a bi-square weighting function so that W(s) is the diagonal matrix with elements

$$\begin{array}{}\text{(3)}& w_{is}={\left(\mathrm{1}-{\left({\displaystyle \frac{d_{is}}{h}}\right)}^{\mathrm{2}}\right)}^{\mathrm{2}}\mathrm{if}{d}_{is}<h\mathrm{and}\mathrm{0}\mathrm{if}{d}_{is}>h.\end{array}$$

Here, h is called the bandwidth. We use the function gwr.basic of the package gw.model on R (Lu et al., 2014). We have chosen h=8 (in longitude and latitude degree units).

The emulator is trained on an ensemble of simulations performed on data sufficiently browsing the plausible input space. The more the input space is filled, the more robust the emulator will be. The fact that we have extended the range of forcing parameters to past, present, and future scenarios is an excellent way to fill this input space. Even if they are not a fully exhaustive sampling, it is reasonable to assume that they represent a model population (Chandler, 2013). The emulator is used to estimate the conditional probability distribution of the ecosystem variables (Table 1) in each location or biome based on the forcing variables (Table 4).

Table 4 The nine global forcing variables and their values for the different periods considered. The global forcing variables are the predictors of the emulator; we give their values for the time slices used in the calibration: GHG atmospheric concentration (carbon dioxide, methane, nitrogen protoxide) in parts per million (ppm), Earth orbit parameters (eccentricity Ecc, obliquity Obl, omega ω), population in millions of people (M), volcanic activity (Volc), and solar activity (Sol, in MeV). Every column is one out of the v forcings of Eq. (1) (columns of X). For columns CO₂, CH₄, N₂O, and population, we use observations until 2010 (see Sect. 2.1) and the values used by the three scenarios (RCP2.6, RCP4.5, RCP8.5) for projections (2020 to 2100).

The total number of simulations available for 18 time slices and three scenarios (for the 2020 to 2100 time slices) from 2 to 18 models is 582. As each one simulates climate for 65 559 grid points of the globe (after downscaling), the total number of available observations is $\sim \mathrm{4}\times {\mathrm{10}}^{\mathrm{7}}$. Because of their strong spatial correlation, it is not necessary to use all the grid points. To have a balance between the different biomes represented on the Earth, for each time slice and each model, we randomly draw a subset of grid points in each biome proportional to its representativity (the proportion of grid points with that biome). So, the calibration is done on about 3×10⁵ observations (the number is slightly lower than expected because some biomes are absent in some simulations).

As the global forcing variables have the same value for all the grid points for a given time slice and scenario, this may produce problems for the computation of the emulator; we have then added to them a small noise value (a Gaussian value with mean 0 and standard deviation equal to the initial standard deviation divided by 100).

Because of collinearity between the predictors (the global forcing variables), their dimensionality is reduced using principal components (PCs) of the normalized variables. The nine variables are reduced into five PCs together explaining 89 % of their total variance. The first PC is mainly related to the greenhouse gases and the global population, which is strongly correlated with them. PC3 shows the opposition between solar and volcanic activity. The other ones are difficult to interpret.

The 13 PFTs are aggregated into eight PFTs, representing the main types of vegetation across the world (Table 5). This enables reducing the dimension of the ecosystem variable vector. The number of variables to be estimated by the emulator is finally reduced to 17 (predictands in Table 6).

Table 5 Aggregation of the plant functional types (PFTs). The 13 original PFTs are defined in Table 1.

Table 6 The t values of the linear regression between the PFT and bioclimatic variables as well as the five PCs of the forcing variables with the R² (in %) of the linear model and the GWR model. For the latter, the coefficients cannot be displayed in a table and are represented as maps in Fig. 5.

For comparison, we first apply a simple linear model to the data. Table 6 shows that most of the bioclimatic variables are not well-fitted. Indeed, only 9 of 17 have a proportion of reconstructed variance higher than 10 % (column R² (lm)), which clearly shows that a single linear model is not adapted to our objective. The GWR model is much better as half of the bioclimatic variables have R² higher than 0.50, but at the cost of a much larger effective number of parameters (ENPs) (ENP ∼ 1311 vs. 6 for the global linear model). As we have a very high number of observations ($\sim \mathrm{3}\times {\mathrm{10}}^{\mathrm{5}}$), the GWRs remain very significant. This is illustrated in Fig. 5 by the maps of the spatial distribution of the regression coefficients (they are standardized by their standard error, which comes to the t value). The degree of smoothing is defined by the bandwidth h (here set to 8). With a lower h, the patterns should be much patchier and the ENPs should become much larger, thereby diminishing the prediction capability of the model.

Figure 5 Local (longitude × latitude) response of a few (bio)climatic variables to the PCs of forcing variables. The colors represent the t values (t-val) of the response (regression coefficient $/$ standard error). Maps 1 to 5 represent the first five PCs. Map 6 represents the R² of the GWRs. Maps 7 and 8 represent the estimates for low-CO₂ time slices (<290 ppm, i.e., LGM and mid-Holocene) and high-CO₂ time slices (>600 ppm, second half of 21st century with scenario RCP8.5), respectively. The last graphic represents the estimates as a function of the observations: blue dots are for observations with CO₂ < 250 ppm, green dots for CO₂ belonging to [250, 400] ppm, orange dots for [400, 600] ppm, and red dots for CO₂ > 600 ppm. The first series of maps is for the net primary production of the ecosystem (NPP_tot); the second is for annual temperature (T_ann), and the last one is for annual precipitation (P_ann).

The interpretation is not straightforward because the predictors are not the forcing parameters but their PCs. Nevertheless, a partial interpretation is possible by using the correlation of the PCs with the forcings. For example, PC1 is correlated with the greenhouse gases (GHGs) and the population size, and Fig. 5 (Tann/PC1) shows that the maximum impact of the GHGs on temperature is in Russia. Figure 5 (Pann/PC1) shows that the GHGs have a negative effect on the precipitation around the Mediterranean and northern Africa and a positive effect on Russia. The decrease in precipitation according to the increase in GHGs is a well-known feature in the Mediterranean (Giorgi and Lionello, 2008). Figure 5 also shows that the annual temperature (Tann/Pred vs. Obs) is well-emulated with a strong correlation with CO₂. It is also true, but to a lesser extent, for NPPtot. The annual precipitation estimates have a much lower variance than the observations.

2.8 Paleodata assimilation (step 4)

The climatic data estimated by the emulator have biases and uncertainties likely larger than those of the ESMs which have generated the calibration data. The first source of errors of the emulator is the same as the source of errors of the ESMs, which are imperfect representations of the real world. The second source of errors comes from the oversimplification, which consists of replacing a complex model by a statistical model, in considering that a few global forcings are able to generate the complex observed climate field. It is partly compensated for by the fact that the emulator is based on a large diversity of model simulations, in contrast to a single ESM. The third source of errors is the variance of the error terms e_j(s) of Eq. (1): the statistical model does not perfectly fit the data. This prompted us to converge the emulator estimates to the real world represented by the paleoclimate data. This is known as paleodata assimilation (Goosse et al., 2012). This led to adjusting a few parameters known to have a large uncertainty to optimize temperature and precipitation using information given by proxies (Table 7).

Table 7 Proxy-based reconstruction of annual temperature (TANN, ^∘C) and annual precipitation (PANN, mm). The values are expressed as anomalies from the preindustrial period for the 10 spatial boxes and the nine time slices, obtained from pollen (Guiot and Kaniewski, 2015) and corrected or made precise as indicated in Table 3.

These uncertain parameters are (i) local impact of volcanic activity, (ii) local impact of solar activity, (iii) a possible global bias for the temperature simulation (δT) and (iv) for the precipitation simulation (δP), and (v) the standard error of temperature (σT) and (vi) of precipitation (σP). Parameters (iii) and (iv) are related to a possible global ESM bias and are then independent of the geographical position. Parameters (v) and (vi) express the mean quadratic difference between observations and simulations and then contain the ESM mean error and the emulator mean error (after removing the effects of biases δT or δP).

The reconstructed annual temperature and precipitation values are expressed as anomalies relative to the present day for each time slice and each box given in Table 7. It shows that the first two time slices (4200 and 3200 yr BP) were characterized by dry conditions in the eastern Mediterranean Basin (boxes Lev, Ana). All the other time slices are characterized by temperature changes.

The statistical simplification of the emulator makes it possible to run it thousands of times in a relatively short computational time as required by the assimilation methods (Widmann et al., 2010). We use a Bayesian approach called the Markov chain Monte Carlo (MCMC) method, which makes it possible to converge towards the best parameters in the sense of probability distribution (Hargreaves and Annan, 2002).

Starting from the prior probability distributions of the six parameters, available climate reconstructions and model outputs provide estimates of their posterior probability distributions. For each time period, these six parameters were optimized to provide the best posterior probability distribution of mean temperature and precipitation at the centers of the 10 Mediterranean boxes. The prior probability distributions are given by uniform distributions delimited by the ±d (Table 8) for the volcanic and solar activity. The prior distribution of the other parameters was assumed to be uniform within wide ranges so that they were considered noninformative for the purpose of this study.

Table 8 Forcing variables used for the 15 time slices and scenarios. The forcing greenhouse gas atmospheric concentrations (CO₂, CH₄, N₂O), orbital parameters (ecc: eccentricity, Obl: obliquity, ω), population size (POP) in number of people, volcanic activity with uncertainty (d), and solar activity with uncertainty (d).

From Tables 3 and 7, we consider the 4200 and 3200 yr BP time slices to provide robust paleoclimatic information for precipitation and the other periods to provide robust information for temperature. Both the precipitation and temperature signals were robust for the present time slice. Even if we ultimately calculate both temperature and precipitation, only the robust variable(s) is (are) used to constrain the emulator.

2.9 The viticulture index (VI)

We subsequently applied our emulator to the question of how viticulture has evolved in the Mediterranean region and in response to which climatic stimuli and global forcing. Numerous bioclimate indices have been published to delimit viticultural zones in the world (Tonietto and Carbonneau, 2004; Santos et al., 2012; Howell, 2001). Among them, the following are cited: (1) the sum of degree days above 10 ^∘C during the growing season or the heliothermal index of Huglin (HI), (2) the number of days with a minimum temperature below −17 ^∘C, which is very important for the grapes growing in continental climates, (3) the minimum temperature of September (cool night index) important for the ripening, (4) the sum of the growing season of the monthly temperature multiplied by monthly precipitation (Hyl), and (5) the drought stress index (DI) related to the potential water balance of the soil during the growing season. Malheiro et al. (2010) proposed a composite index (CompI) calculated based on the ratio of years simultaneously satisfying four criteria (HI > 1400 degree days), DI > −100 mm, Hyl < 5100 ^∘C mm⁻¹, and Tmin > −17 ^∘C).

Some of the climate variables needed for these indices were not available from the BIOME4 outputs. However, other variables, such as those associated with the net primary production of plant types, are very interesting because they include the CO₂ effect on photosynthesis and water use efficiency. Considering only rain-fed viticulture, we propose the following index, VI:

$$\begin{array}{}\text{(4)}& \mathrm{VI}=\left(\mathrm{1}-I_{\mathrm{NPPtrop}}\right)\cdot I_{\mathrm{NPP}}\cdot I_{\mathrm{Pann}}\cdot I_{\mathrm{MTWA}}\cdot I_{\mathrm{MTCO}}\cdot I_{\mathit{\alpha }},\end{array}$$

where each of these factors denoted I_x follows the function

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{I}_x=\mathrm{0}\mathrm{for}x<x_{\min }\\ {\displaystyle }& {\displaystyle }{I}_x=\mathrm{1}\mathrm{for}x>x_{\max }\\ {\displaystyle }& {\displaystyle }{I}_x={\displaystyle \frac{x-x_{\min }}{\left(x_{\max }-x_{\min }\right)}}\mathrm{for}{x}_{\min }\le x\le x_{\max }.\end{array}$$

For x=NPP_trop, the net primary production of tropical plants, the interval [x_min, x_max] is [0, 10 kg C m⁻²]. For x=NPP, the total ecosystem net primary production, it is [500, 1000 kg C m⁻²]. For x=P_ann, the total annual precipitation, it is [400, 800 mm]. For x=MTWA, the mean temperature of the warmest month, it is [18, 23 ^∘C]. For x=MTCO, the mean temperature of the coldest month, it is [3, 12 ^∘C]. The lower limit of 3 ^∘C corresponds to a minimum mean monthly temperature with no chance of having an absolute daily minimum below −17 ^∘C. For x=α, the actual to potential evapotranspiration ratio, it is [30 %, 60 %].

All factors except those related to NPP_trop assume that vine growth is limited only by the lower values and not by the higher values. Because there is no possible viticulture in the tropics (where temperature is not cold enough for an appropriate dormancy), $\mathrm{1}-I_{{\mathrm{NPP}}_{\mathrm{trop}}}$ is limited by its upper value.

This approach is a static approach as it is based on an equilibrium vegetation model only based on the mean state of the climate. The variability of the climate and therefore the extreme values of some variables are not considered. We are aware of these limitations, and it is the reason why our interpretations will only concern the potential extension of viticulture. Other weaknesses of the approach are the impossibility to consider the diversity of cultivars, the adaptation of grapes to frost, and the risk of biotic damage to the grapevines. This is a niche model approach of the same nature as those used in the literature cited at the beginning of this section. This oversimplified approach nevertheless has the advantage of providing a framework for the past and future trends and of being easily applicable to the past for which the data are rather limited.

3 Results

3.1 Paleodata assimilation

From the posterior distributions (Table 9), it appears that the simulated volcanic activity impacts were the strongest during the cold periods (2500, 1300, 700, and 200 yr BP) and rather weak during the dry periods (4200 and 3200 yr BP) and warm periods (2000 and 1000 yr BP). The confidence intervals are significant as they do not overlap. For the present, volcanic and solar activities have no significant impact, which is easily explained by the dominant GHG effect (not used for constraining the assimilation). The impact of solar activity is not clear, as there is no significant difference between cold and warm periods or between the dry and wet periods. The temperature bias (δT independent of the spatial variability) is estimated to be between 0.6 and 1.6 ^∘C for the periods between 2500 and 200 yr BP and is not significantly different from zero for the dry periods and for the present. The δP estimates were not significant for any time period. Figure 6 presents the overall correlations between the emulator outputs and the proxy-based reconstructions. A good emulation should be indicated by good coherency of both blue and black lines and a low variance of the dots around the black line. So, the temperature was particularly well-simulated by the emulator with a squared correlation (R²) of 0.76, but with an overestimation of the low temperatures. As expected, precipitation is less well-simulated with a significant R² of 0.22 with an underestimation of the large (negative or positive) anomalies.

Table 9 Posterior distribution of the parameters. The estimates are given by the best fit, and the CI is the 90 % confidence interval calculated on the 10 % best fits. Volc index has no units; the solar index is in megaelectron volts (MeV), T is in degrees Celsius (^∘C), and P is in millimeters per year (mm yr⁻¹).

Figure 6 Temperature and precipitation anomalies for the 10 Mediterranean boxes estimated by data assimilation versus data reconstructed from proxies. Temperature dots represent the mean annual temperature in the 10 boxes for the 11 periods from 2500 yr BP until the present, and precipitation dots represent the total annual precipitation in the two oldest periods (4200 and 3200 yr BP) and in the present. The vertical bars are the 90 % confidence intervals (CIs). The blue line is the perfect estimate line (R=1), and the black line is the regression line between the simulations and proxy reconstructions.

Figure 7 Spatial distribution of the temperature and precipitation anomalies for the past and present time slices. The pairs of symbols are simulations (left symbols) and observations or pollen reconstructions (right symbols). The colors are used for the temperature anomalies. Triangles are used for negative precipitation anomalies and circles for the positive precipitation anomalies. The size of the triangles and circles is proportional to the absolute value of the precipitation anomalies.

3.2 Application of the emulator to past and future scenarios

This assimilation process was applied to simulate the annual temperature and precipitation for the 10 boxes and the 13 time slices and scenarios. The spatial patterns of the reconstructions and simulations are consistent, except for the temperatures at 4200 and 3200 yr BP. (Fig. 7). Therefore, global forcing variables can drive not only the mean climate evolution but also the spatial patterns reconstructed by the proxy data in the Mediterranean. Volcanic activity seems to be the main driver for the cold periods of the Iron Age (2500 yr BP), the LALIA (1300 yr BP), the LIA (700, 200 yr BP), and the warm periods of the RCO (2000 yr BP) and MCA (1000 yr BP). The main driver of the present warming and drying (present time in Fig. 7) is GHG forcing. For the 4200 and 3200 yr BP periods, the volcanic and solar activity drivers do not seem to be plausible explanations for the droughts, at least according to our emulator, and likely the ESMs; this also explains the absence of agreement.

Future projections (Fig. 8) show general warming for the four scenarios (+1.5C, +2C, +3C, +5C) that is particularly strong for +3C and +5C scenarios. The local warming is less strong in the areas influenced by the Atlantic Ocean (France and the Iberian Peninsula). For precipitation, the signal is weak for +1.5C and +2C with drier and wetter zones and clearer for +3C and +5C (all the symbols are triangles, i.e., negative). Combined with the warming, it is undeniable that water stress will increase considerably with +3C and +5C scenarios.

Figure 8 Spatial distribution of the temperature and precipitation anomalies simulated for future time slices (scenarios). The colors are used for the temperature anomalies. Triangles are used for negative precipitation anomalies and circles for positive precipitation anomalies. The size of the triangles and circles is proportional to the absolute value of the precipitation anomalies.

To answer the question of whether volcanic activity can mitigate the impact of high GHG emissions, Fig. 8 shows that the temperature distributions of +5 CV+ and 5 CV− are quite similar, meaning that the cooling effect of volcanoes is relatively low in a large GHG emission scenario. The precipitation distributions of +5 CV+ and 5 CV− are similar in the western Mediterranean, and +5 CV+ produced slightly higher precipitation anomalies than +5 CV− in the eastern Mediterranean.

3.3 Independent validation

Additional independent validation was completed through comparisons with a tree-ring-based reconstruction of the Palmer Drought Severity Index (PDSI) (Cook et al., 2015) for the last 2 millennia, for which tree-ring data are available. The PDSI reflects the spring–summer soil moisture conditions. We compared this variable with the reconstruction of the $E/\mathrm{PE}$ (ratio of actual to potential evapotranspiration in %) variable provided by BIOME4. $E/\mathrm{PE}$ is a moisture index which is equal to zero when the soil is fully dry and 100 when it is fully wet. The range of the PDSI is usually between −6 and 6 units. Negative values correspond to conditions drier than normal. Considering that both indices are not fully similar, a visual comparison shows pretty good agreement (Fig. 9). So for 2000 BP, Iberia and central Europe were wet in both maps, and no PDSI data were available for the southern and eastern Mediterranean. For 1300 BP, all of Europe was dry in both maps except eastern Spain, which was as wet as at 2000 BP. For 1000 BP, the area surrounding the sea was wet except Greece in both maps, and for 700 BP all areas were wetter than 1000 BP in both maps.

Figure 9 Comparison of $E/\mathrm{PE}$ reconstructed in this paper with independently reconstructed Palmer Drought Severity Index (PDSI). The PDSI anomalies are based on the 1928–1978 reference period and are reconstructed from tree rings (Cook et al., 2015). $E/\mathrm{PE}$ is the ratio of actual to potential evapotranspiration obtained from the emulator (in %).

3.4 Evolution of viticulture from the Bronze Age to the end of the 21st century

Applied to the mean climate of 1980–2009, the viticulture index (VI) approximately reproduces the area in Europe where viticulture is present (Fraga et al., 2013) (Fig.10). As shown in Fig. 10a, the simulations give approximately the same limits, with some noticeable exceptions in Germany as well as the Alpine and Balkan regions, which are likely too continental to pass the MTCO criterion. This is the main limitation of our analysis. The presence of viticulture greatly depends on local factors such as valleys, slopes, and soil, which are not accounted for in our analysis. More generally, these discrepancies show that our index is not calibrated for continental regions with winters colder than 3 ^∘C (on average). This threshold of 3 ^∘C is likely overestimated by about 5 ^∘C. So, we consider this index to adequately represent the macroclimate of potential vine distribution in western Europe and the Mediterranean region.

Figure 10 Distribution of the wine sites: (a) European wine regions (circles) with the corresponding CompI value for the period of 1980–2009: from Fig. 1e of Fraga et al. (2013). (b) Simulation using the viticulture index calculated using 1961–1990 climate data.

The variables needed for viticulture in Eq. (4) are provided by the emulator for all time slices studied in the past time slices and for future scenarios. The procedure is sketched in Fig. 11. We modified the global forcings according to the values of Table 3 for the past and Table 2 for the future, and we thereby simulated, by applying the emulator, viticulture extension maps (Fig. 12). Figure 12 shows that during the dry time slices of 4200 and 3200 yr BP, the suitable areas were located between 34 and 46^∘ N latitudes. During the cold periods (2500, 1300, 700, and 200 yr BP), they occupied approximately the same area between 34 and 47^∘ N. During the warmer periods (2000 and 1000 yr BP), the southern limit did not change much, but the northern limit reached 49^∘ N, implying that most of Gaul was suitable for viticulture, as already shown by Bernigaud et al. (2021). The present map is warmer than the preindustrial period (200 yr BP slice), which suggests that viticulture is now at its maximum potential extension in 4000 years, up to 51^∘ N. Because of the drier conditions, the southern limit has shifted from 34 to 35^∘ N. Note that these variations only depend on climate changes, without any consideration of other factors such as types of soils, which are very important for wine quality.

Figure 11 Diagram of the application of the emulator to reconstruct and predict the viticulture extension.

Figure 12 Distribution of the viticulture index (Vit. Index) for several time slices of the past (in yr BP) and future scenarios (expressed as anomalies of global temperature vs. preindustrial reference). The index is labeled VI and is explained in Eq. (4). White areas indicate where it is impossible to cultivate vine.

In the future projections, the northern extension of potential viticulture should shift from 51^∘ N (+1.5C scenario) to 53^∘ N (+3C scenario) and even above 55^∘ N (for the +5C scenario). This would allow viticulture to be possible up to central England, but it would regress in the south due to stronger water stress. In the Iberian Peninsula, only the Atlantic coast should be suitable for cultivating wine grapes unless there is significant irrigation. These projected unfavorable conditions were in agreement with Fraga et al. (2013) based on other viticulture indices. The 5 CV+ and 5 CV− scenarios appear quite like the 5C scenario, indicating that the effect of high or low volcanic activities should have a weak effect on the potential distribution of viticulture in comparison to a strong GHG forcing.

The results are summarized in Fig. 13. In the past, only the warm regions of the southern band (29–37^∘ N) had a suitable wine-growing area equivalent to the present. The central geographical band (37–44^∘ N) and the northern one (44–48^∘ N) underwent sudden changes in the viticulture area, first from the cool Iron Age (2500 yr BP) to the RCO (2000 yr BP) and later from the end of the LIA (200 yr BP) to the present. The cold periods are all characterized by a decline in viticulture at latitudes above 37^∘ N. For the future projections, northward displacements are likely to be drastic from +3 ^∘C global warming. Viticulture potential will likely disappear from North Africa and is set to decrease drastically in the Iberian Peninsula. In contrast, potential productive areas will likely expand at latitudes above 50^∘ N, in the Balkans, and in the Alps, which is a significant change considering that our approach underestimates the grapevine potential of continental climates.

Figure 13 Synthesis of the evolution of viticulture in the Mediterranean area. The curves represent three bands of latitude: the southern band with latitude < 37^∘ N in red, the center band with latitudes between 37 and 44^∘ N in green, and the northern band with latitudes between 44 and 48^∘ N in blue). The thin lines show the results for the boxes corresponding to the Iberian Peninsula, Anatolia, and France (definition in Fig. 4). The scenarios and time slices are defined in Table 8. The x axis gives the center of the time slices considered or the future scenario.

For the two additional scenarios of high emissions combined with extreme volcanic and solar activities, the question is whether an entirely hypothetical set of strong volcanic events could slow down the decline in viticulture in the south. The answer is slightly positive for Spain, Italy, Greece, and Turkey (green curve increases for +5 CV+ and decreases for +5 CV−), but it is negative in North Africa and the Levant (red curve) because the water stress should remain too substantial. The effect is also positive in the northern band and Turkey (the blue and cyan curves increase for +5 CV+ and decrease for +5 CV−). In all the cases, the effect was clearly too small to compensate for the GHG effect.

4 Discussion

The discussion considers six points, which concern methodological innovations as well as contributions to climate science and viticulture science.

1. We used an emulator approach which recently became popular in the ESM community. In some ways, it has also been used with paleoclimatic data (Crucifix, 2012; Joos et al., 1996),= and to emulate a climate–vegetation system (Foley et al., 2016). But this study is the first to use an ensemble of climate models, along with several forcing variables provided by past, present, and future states, to project vegetation conditions from global climate drivers. This has involved solving several technical issues.

2. We showed that the major climatic variations of the last millennia in the Mediterranean Basin can be attributed to volcanic activity, whereas the effects of solar activity were negligible. The effects of volcanic and solar activities have been largely debated, as the reconstruction of temperatures in the last millennium from tree rings has often shown less significant cooling than the model simulations (Mann et al., 2012). A more recent tree-based reconstruction (Stoffel et al., 2015) showed substantial summer cooling after the Samalas eruption in 1257 and that of Tambora in 1816, but less intense than simulated. Luterbacher et al. (2016) showed that solar activity had a relatively weak influence on the European summer temperatures. Thus, our results are in line with the state of the art in the scientific literature.

3. The effect of this volcanic forcing has a clear spatial pattern across the Mediterranean Basin. From 2500 yr BP to the present, temperature variations were more significant in the north and in the west than in the southeast (Fig. 7). Fischer et al. (2007) found that northern and western Europe were the coldest and driest areas after an eruption. They also found that southern Europe, North Africa, and the Levant have experienced milder and wetter weather than at present. Part of this pattern might be due to regional feedbacks of vegetation on the climate. Some climate model simulations have shown (Reale and Dirmeyer, 2000) that wetter vegetation during the RCO may have induced a northward shift of the Intertropical Convergence Zone during the summer over the African continent with an increase in moisture in North Africa and Iberia, as well as a decrease in the central Mediterranean (Reale and Shukla, 2000). In the future, the southeast should be relatively less dry and warmer than the northwest, especially in summer, which is consistent with our results (Fig. 8) (Giorgi and Lionello, 2008).

4. Our simulations of climate change effect on viticulture are mostly concerned with larger-scale regional production systems since they would require nearly full-time engagement of wine growers with high certainty of production every year, sufficient for speculative trade. For smaller domains and local consumption, it may have also been possible to produce wine under degraded or unstable weather conditions. For example, in England, viticulture continued to be practiced on land owned by the church, even as risks increased due to a wetter climate, with cooler summers and milder winters (Clout, 2013). In most wine regions in western Europe, particularly in France, the grape harvest dates were advanced after the LIA and particularly after the 1940s (Le Roy Ladurie et al., 2006). For example, from 1945 to the beginning of the 21st century, in Chateauneuf du Pape (southern Rhone Valley, France) the harvest date advanced on average from 1 October to 11 September (Ganichot, 2002). This change is related to summer warming, but factors related to wine quality, agricultural practices, and alcohol content regulation may induce a bias in the interpretation (de Cortázar-Atauri et al., 2010). In the Languedoc region, the potential alcohol content increased by two degrees from 1984 to 2013 (van Leeuwen and Darriet, 2016). Even if the alcohol content does not exactly reflect the grape yields, earlier harvest dates with a higher sugar content have clearly been related to improved conditions of grape cultivation since the LIA, which is also related to increased productivity. Another climate-sensitive symbol of Mediterranean agriculture is olive trees. Moriondo et al. (2013) found three characteristic periods in olive growing, namely the RCO (300 BCE to 400 CE) and the MCA (900 CE to 1200 CE), during which olive-growing areas expanded northward, and the LIA (1400 CE to 1900 CE) during which a contraction was reported. This is consistent with the variations in viticulture, at least in western Europe and the Mediterranean area as our VI is not properly calibrated for strongly continental climates.

5. Major difficulties are forecast for 21st century viticulture in Spain and North Africa. These are particularly important for global warming levels of +3 ^∘C and more. Caffarra and Eccel (2011) showed that the projected warming should make some mountain sites at approximately 1000 m climatically suitable for viticulture before the end of this century. The MCA limit will certainly be passed. However, other factors could become limiting, such as excessively mild winters that enable pest attacks and infections, lack of expertise in vine growing and wine making, and products that are more expensive than the current Mediterranean wines (Clout, 2013). Other limitations are extreme events (late frost, flooding). More frequent and more intense heatwaves will no longer be favorable to viticulture at the present southern Mediterranean limit of its niche. These factors are not considered in our approach. A recent study considering grape varieties showed (Sgubin et al., 2022) that, for global temperature anomalies below 2 ^∘C, the mean relative area loss was estimated to be 3.9 % ^∘C⁻¹, while for higher values of global warming this loss trend is estimated to be a much larger rate of 17.1 % ^∘C⁻¹. This confirms the existence of a safe limit below 2 ^∘C of global warming for the European wine-making sector, above which adaptation might become far more challenging.

6. It is not very likely that intense volcanic activity could slow down this decline because (despite the many other negative consequences of such events) the beneficial climatic effects of this activity would be highest in regions with increased wine-growing potential under global warming (Turkey, northern Europe, the Alps, and the Balkans) and negligible in North Africa. IPCC has assessed the literature concerning the question of whether volcanic eruptions could be analogous for geoengineering proposals for climate mitigation (Myhre et al., 2013). Independent of their possible negative side effects, solar radiation modification (SRM), consisting of injecting sulfur aerosols in the stratosphere like a volcano does, is likely not an efficient technique to cool the atmosphere. Indeed, our results show that very strong volcanic activity, even accompanied by low solar activity, should not have any significant effects on viticulture extension in comparison to the radiative forcing from anthropogenic GHG emissions.

5 Conclusion

We demonstrated the efficiency of a statistical emulator based on multiple ESMs and calibrated over a large range of forcing and climate states to link the production of a regionally important crop (grapes) to global climate forcing variables. There are two main methodological innovations: (i) past climate simulations are used together with future ones to calibrate a robust emulator, and (ii) this emulator is independent of any single model because it captures what is common among all the available model simulations. But it remains an emulator which captures no more than the information contained in the ESMs.

Using it on past time slices, we showed that volcanic activity is a good predictor of the past temperature variations in the Mediterranean Basin and consequently of the viticulture productivity. During the warm phases of historical times (the Roman Climate Optimum and the Medieval Climate Anomaly) characterized by low volcanic activity, the viticulture area has shifted northward from 47 to 49^∘ N. This historical limit has already been passed at the present time as it has already shifted to 51^∘ N, and with global warming of +3 ^∘C, it should pass the 53^∘ N limit. Even if, in the warm periods of the past, North African and Iberian farmers never had real difficulties in cultivating grapes, they will meet large difficulties with global warming of +3 ^∘C or more, except on the Atlantic margin. Moreover, our sensitivity experiments show that even intense volcanic activity is not sufficient to stop this decline.