# Using climate change scenarios to simulate mobility of metal contaminants in soils: the example of copper on a European scale

Laura Sereni, Julie-Mai Paris, Isabelle Lamy, Bertrand Guenet

► To cite this version:

Laura Sereni, Julie-Mai Paris, Isabelle Lamy, Bertrand Guenet. Using climate change scenarios to simulate mobility of metal contaminants in soils: the example of copper on a European scale. Soil, 2023, 2023, pp.2350. 10.5194/egusphere-2023-2350. hal-04535489

HAL Id: hal-04535489

https://hal.science/hal-04535489

Submitted on 6 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

1 Using climate change scenarios to simulate mobility of metal contaminants in soils: the example of copper on

2 a European scale

3

4 Laura SERENI¹,², Julie-Maï PARIS³, Isabelle LAMY¹, Bertrand GUENET³

5 ¹Université Paris-Saclay, INRAE, AgroParisTech, UMR EcoSys, 91120 Palaiseau, France

6 ²Present address: Univ. Grenoble Alpes, CNRS, INRAE, IRD, Grenoble INP, IGE, Grenoble, France

7 ³Laboratoire de Géologie ENS, PSL Research University, CNRS, UMR 8538, IPSL, Paris, France

8 *Correspondence to Laura Sereni (laurasereni@yahoo.fr)

9

10 Abstract:

11 Soil contaminant deposition is highly dependent on anthropogenic activities while contaminant retention,
12 mobility and availability are highly dependent on soil properties. The knowledge of partitioning between soil solid
13 and solution phases is necessary to estimate whether deposited amounts of contaminants will rather be leached
14 through runoff or accumulated. Besides pedological driven partitioning, runoff is expected to change during the
15 next century due to changes in climate and in rainfall patterns. In this study, we aimed at estimating at the
16 European scale the areas concerned by potential risk due to contaminant leaching (LP). We also defined in the
17 same way the surface areas where limited Cu leaching occurred, leading to potential accumulation (AP) areas.
18 Among contaminant, we focused on copper (Cu) widely used in agriculture, resulting in high spatial variations in
19 deposited and incorporated amounts in soils. We developed a method using both Cu partition coefficients (Kf)
20 between total and dissolved Cu forms, and runoff simulation results for historical and future climates. The
21 calculation of Kf with pedo-transfer functions allowed us to avoid any uncertainties due to past management or
22 future depositions that may affect total Cu concentrations. Areas with high potential risk of leaching or of
<figure><img src="image_2.png" /></figure>

1 Using climate change scenarios to simulate mobility of metal contaminants in soils: the example of copper on

2 a European scale

3

4 Laura SERENI¹,², Julie-Maï PARIS³, Isabelle LAMY¹, Bertrand GUENET³

5 ¹Université Paris-Saclay, INRAE, AgroParisTech, UMR EcoSys, 91120 Palaiseau, France

6 ²Present address: Univ. Grenoble Alpes, CNRS, INRAE, IRD, Grenoble INP, IGE, Grenoble, France

7 ³Laboratoire de Géologie ENS, PSL Research University, CNRS, UMR 8538, IPSL, Paris, France

8 *Correspondence to Laura Sereni (laurasereni@yahoo.fr)

9

10 Abstract:

11 Soil contaminant deposition is highly dependent on anthropogenic activities while contaminant retention,
12 mobility and availability are highly dependent on soil properties. The knowledge of partitioning between soil solid
13 and solution phases is necessary to estimate whether deposited amounts of contaminants will rather be leached
14 through runoff or accumulated. Besides pedological driven partitioning, runoff is expected to change during the
15 next century due to changes in climate and in rainfall patterns. In this study, we aimed at estimating at the
16 European scale the areas concerned by potential risk due to contaminant leaching (LP). We also defined in the
17 same way the surface areas where limited Cu leaching occurred, leading to potential accumulation (AP) areas.
18 Among contaminant, we focused on copper (Cu) widely used in agriculture, resulting in high spatial variations in
19 deposited and incorporated amounts in soils. We developed a method using both Cu partition coefficients (Kf)
20 between total and dissolved Cu forms, and runoff simulation results for historical and future climates. The
21 calculation of Kf with pedo-transfer functions allowed us to avoid any uncertainties due to past management or
22 future depositions that may affect total Cu concentrations. Areas with high potential risk of leaching or of
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

accumulation were estimated over the XXIth century by comparing Kf and runoff to their respective European
median. Thus, at three distinct times, we considered a grid point at risk of LP if its Kf was low compared to the
European median and its runoff was high compared to the European median of the time. Similarly, a grid point
was considered at risk of AP if its Kf was high and its runoff was low compared to their respective European median
of the time. To deal with uncertainties in climate change scenarios and the associated model projections, we
performed our study with two representative atmospheric greenhouse gases concentration pathways, defined
with climate change associated to a large set of socio-economic scenarios found in the literature. We used two
land surface models (ORCHIDEE and LPJmL, given soil hydrologic properties) and two global circulation models
(ESM2m and CM5a, given rainfall forecast). Our results show that, for historical scenario 6.4 ± 0.1 % (median,
median deviation) and 6.7 ± 1.1 % of the grid cells of the European land surfaces are concerned by LP and AP
respectively. Interestingly, our results simulate a constant global surface concerned by LP and AP, around 13% of
the grid cells, consistent with an increase in AP and a decrease in LP. Despite large variations in LP and AP extents
depending on the land surface model used for estimations, the two trends were more pronounced with RCP 6.0
than with RCP 2.6, highlighting the global risk of combined climate change and contamination and the need for
more local assessment. Results are discussed to highlight the points requiring improvement to refine predictions.

Keywords: regional modeling, transfer functions, ISIMIP, LUCAS Topsoil data, mapping risk
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

accumulation were estimated over the XXIth century by comparing $K_f$ and runoff to their respective European
median. Thus, at three distinct times, we considered a grid point at risk of LP if its $K_f$ was low compared to the
European median and its runoff was high compared to the European median of the time. Similarly, a grid point
was considered at risk of AP if its $K_f$ was high and its runoff was low compared to their respective European median
of the time. To deal with uncertainties in climate change scenarios and the associated model projections, we
performed our study with two representative atmospheric greenhouse gases concentration pathways, defined
with climate change associated to a large set of socio-economic scenarios found in the literature. We used two
land surface models (ORCHIDEE and LPJmL, given soil hydrologic properties) and two global circulation models
(ESM2m and CM5a, given rainfall forecast). Our results show that, for historical scenario 6.4 ± 0.1 % (median,
median deviation) and 6.7 ± 1.1 % of the grid cells of the European land surfaces are concerned by LP and AP
respectively. Interestingly, our results simulate a constant global surface concerned by LP and AP, around 13% of
the grid cells, consistent with an increase in AP and a decrease in LP. Despite large variations in LP and AP extents
depending on the land surface model used for estimations, the two trends were more pronounced with RCP 6.0
than with RCP 2.6, highlighting the global risk of combined climate change and contamination and the need for
more local assessment. Results are discussed to highlight the points requiring improvement to refine predictions.

Keywords: regional modeling, transfer functions, ISIMIP, LUCAS Topsoil data, mapping risk
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

39

40

41
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

39

40

41
<figure><img src="image_6.png" /></figure>

42 1. Introduction

43

44 At a large spatial scale, trace element contents in soils are highly variable in relation with the trace element
45 contents of the soil parental rocks and with local anthropogenic inputs of various origins (Flemming and Trevors,
46 1989; Noll, 2003; Salminen and Gregorauskiene, 2000). Some trace elements like copper (Cu) or zinc are required
47 for several biological mechanisms, but when highly concentrated they may have toxic effects on soil organisms
48 (Giller et al., 1998). In particular, Cu is widely used, as a fungicide, especially against downy mildew in vineyards
49 (Komárek et al., 2010), but also in industrial processes. At the European scale, a gradient of soil Cu concentrations
50 can be find from typical baseline values between 5 mgCu.kg⁻¹ to 20 mgCu.kg⁻¹ (Salminen and Gregorauskiene,
51 2000), to values larger than 100 mgCu.kg⁻¹, common in cultivated soils and especially in vineyards parcels (Ballabio
52 et al., 2018). It is commonly accepted to conceptually partition the total soil Cu content into different pools of Cu
53 forms in close equilibrium. Briefly, three pools can be defined: a so-called ‘inert’ pool corresponding to Cu included
54 into minerals, a so-called ‘labile’ pool corresponding to Cu sorbed to soil constituents but that can be mobilized
55 according to environmental conditions, and a smallest ‘mobile’ pool corresponding to Cu in soil solution that may
56 be readily available for living organisms but also for transport within soil horizons (Broos et al., 2007; Rooney et
57 al., 2006; West and Coombs, 1981). Schematically, these pools are governed by processes like exchange,
58 complexation or sorption. Also, local soil characteristics such as organic matter, pH or cationic exchange capacity
59 can affect the proportion of Cu in these different pools (Vidal et al., 2009). Any modifications in soil properties or
60 soil solution composition may thus affect Cu equilibrium between sorbed and solution phases. The pool of Cu in
61 the solution phase can be assimilated to a potential pool of Cu leaching. Conversely, Cu bound to the solid phases
62 can be assimilated to a potential pool of Cu accumulation in soil. Depending on the main process involved, for a
63 given amount of Cu deposited on soil, the proportions of leached and accumulated Cu can varied from place to
64 place and with time. However, studies simulating whether the soil will rather leach or accumulate a contaminant,
65 i.e. will act as a source or a sink for the contamination, are scarce especially at a large spatial scale. This knowledge,
66 however, could allow to highlight contaminated areas with a potential to leach, disperse or accumulate
<figure><img src="image_7.png" /></figure>

42 1. Introduction

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

At a large spatial scale, trace element contents in soils are highly variable in relation with the trace element contents of the soil parental rocks and with local anthropogenic inputs of various origins (Flemming and Trevors, 1989; Noll, 2003; Salminen and Gregorauskiene, 2000). Some trace elements like copper (Cu) or zinc are required for several biological mechanisms, but when highly concentrated they may have toxic effects on soil organisms (Giller et al., 1998). In particular, Cu is widely used, as a fungicide, especially against downy mildew in vineyards (Komárek et al., 2010), but also in industrial processes. At the European scale, a gradient of soil Cu concentrations can be found from typical baseline values between 5 mgCu.kg⁻¹ to 20 mgCu.kg⁻¹ (Salminen and Gregorauskiene, 2000), to values larger than 100 mgCu.kg⁻¹, common in cultivated soils and especially in vineyards parcels (Ballabio et al., 2018). It is commonly accepted to conceptually partition the total soil Cu content into different pools of Cu forms in close equilibrium. Briefly, three pools can be defined: a so-called 'inert' pool corresponding to Cu included into minerals, a so-called 'labile' pool corresponding to Cu sorbed to soil constituents but that can be mobilized according to environmental conditions, and a smallest 'mobile' pool corresponding to Cu in soil solution that may be readily available for living organisms but also for transport within soil horizons (Broos et al., 2007; Rooney et al., 2006; West and Coombs, 1981). Schematically, these pools are governed by processes like exchange, complexation or sorption. Also, local soil characteristics such as organic matter, pH or cationic exchange capacity can affect the proportion of Cu in these different pools (Vidal et al., 2009). Any modifications in soil properties or soil solution composition may thus affect Cu equilibrium between sorbed and solution phases. The pool of Cu in the solution phase can be assimilated to a potential pool of Cu leaching. Conversely, Cu bound to the solid phases can be assimilated to a potential pool of Cu accumulation in soil. Depending on the main process involved, for a given amount of Cu deposited on soil, the proportions of leached and accumulated Cu can varied from place to place and with time. However, studies simulating whether the soil will rather leach or accumulate a contaminant, i.e. will act as a source or a sink for the contamination, are scarce especially at a large spatial scale. This knowledge, however, could allow to highlight contaminated areas with a potential to leach, disperse or accumulate
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

67 contaminants, and therefore help for long term environmental management.

68 Concurrently, climate change due to anthropogenic activities is expected to impact rainfall patterns in the
69 forthcoming the decades, leading to changes in the frequency and intensity of weather events at regional and
70 local levels (Christensen and Christensen, 2003). For instance, projections forecast an increase in rain- and snow-
71 fall events in winter in Northern Europe but a decrease in summer in the Mediterranean region, which extends to
72 northward regions (Douville et al., 2021). The extent of rain- and snow-fall alterations depends on anthropogenic
73 activities and associated climate change. Thus, climate change will alter water flows throughout the century
74 (Mimikou et al., 2000). For instance, increase in rainfall intensity and in water accumulation in the soil surface due
75 to limited water infiltration may induce large runoff (Chu et al., 2019). Changes in runoff will also change fluxes
76 of elements or of particulates in the soil solution as it has been shown for Cu (Babcsányi et al., 2016). However,
77 the relationships between these changes in runoff and fluxes of elements is still poorly predicted for the next
78 decades.

79 In this framework, our aim was twofold: i) estimate the soil Cu leaching potential areas in Europe, thereafter,
80 named LP, for the beginning of the century and ii) predict their evolution according to different climate change
81 scenarios. Additionally, we aimed to estimate the Cu accumulation potential areas thereafter named AP. We
82 hypothesized that the role of soil as Cu sink or source, linked to the processes of Cu accumulation or leaching, can
83 be described by the combined effects of local runoff amounts and of local soil properties controlling the partition
84 of total Cu in sorbed and solution species. Due to the lack of information about the future Cu deposition and on
85 soil Cu concentrations whatever its form, we developed a method using the partition coefficient (Kf) at the
86 equilibrium between solid and solution phases to determine areas with high or low potential of leaching whatever
87 total Cu concentration. Regarding the lack of data about future deposited amount at large scale, using Kf was
88 necessary to estimate the Cu mobility potential. The LP or AP areas were thus estimated through the combined
89 use of Kf, calculated with the help of pedo-transfer functions, and the use of soil runoff amounts extracted from
90 earth system simulations. With the use of Kf we avoided the uncertainties due to past land management and
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

67 contaminants, and therefore help for long term environmental management.

68 Concurrently, climate change due to anthropogenic activities is expected to impact rainfall patterns in the
69 forthcoming the decades, leading to changes in the frequency and intensity of weather events at regional and
70 local levels (Christensen and Christensen, 2003). For instance, projections forecast an increase in rain- and snow-
71 fall events in winter in Northern Europe but a decrease in summer in the Mediterranean region, which extends to
72 northward regions (Douville et al., 2021). The extent of rain- and snow-fall alterations depends on anthropogenic
73 activities and associated climate change. Thus, climate change will alter water flows throughout the century
74 (Mimikou et al., 2000). For instance, increase in rainfall intensity and in water accumulation in the soil surface due
75 to limited water infiltration may induce large runoff (Chu et al., 2019). Changes in runoff will also change fluxes
76 of elements or of particulates in the soil solution as it has been shown for Cu (Babcsányi et al., 2016). However,
77 the relationships between these changes in runoff and fluxes of elements is still poorly predicted for the next
78 decades.

79 In this framework, our aim was twofold: i) estimate the soil Cu leaching potential areas in Europe, thereafter,
80 named LP, for the beginning of the century and ii) predict their evolution according to different climate change
81 scenarios. Additionally, we aimed to estimate the Cu accumulation potential areas thereafter named AP. We
82 hypothesized that the role of soil as Cu sink or source, linked to the processes of Cu accumulation or leaching, can
83 be described by the combined effects of local runoff amounts and of local soil properties controlling the partition
84 of total Cu in sorbed and solution species. Due to the lack of information about the future Cu deposition and on
85 soil Cu concentrations whatever its form, we developed a method using the partition coefficient (Kf) at the
86 equilibrium between solid and solution phases to determine areas with high or low potential of leaching whatever
87 total Cu concentration. Regarding the lack of data about future deposited amount at large scale, using Kf was
88 necessary to estimate the Cu mobility potential. The LP or AP areas were thus estimated through the combined
89 use of Kf, calculated with the help of pedo-transfer functions, and the use of soil runoff amounts extracted from
90 earth system simulations. With the use of Kf we avoided the uncertainties due to past land management and
<figure><img src="image_9.png" /></figure>

previous Cu deposition and focused on risks arising from future deposition. To do so, we first reviewed the empirical equations estimating Cu's Kf based on soil properties to highlight generic soil properties governing this partition. From this review, we extracted the best compromise Kf equation to estimate partitioning at the regional scale, which ensures more accurate Kf calculation based on pedo-geochemical data typically recorded in soil surveys, thus mainly available. This allowed us to estimate Cu's Kf values to be used at the European scale based on pedo-geochemical soil surveys without the knowledge of soil Cu total content. We then focused on the current state of the climate and its projected changes over the XXIst century, based on two climate change scenarios. To capture the difficulties in runoff prediction and to disentangle the uncertainties between rainfall prediction and runoff calculations of land surface models, we used a set of simulations provided by the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP). These simulations used different land surface models driven by different climate forcings computed by different climate models. For each scenario and each couple of land surface model and climate forcing we estimated the LP or AP of each grid points by comparison between the local values of Kf and of runoff to the respective calculated European median.

## 2. Materials and methods

### 2.1. Equations to estimate copper Kf

The rigorous definition of Kf is based on the concentration ratio of sorbed vs solution species (here Cu) at the equilibrium. Yet, for practical reasons of measurement and applicability, Kf is conventionally derived from total Cu and not from sorbed Cu (Degryse et al., 2009). A general form of the Cu partition coefficient between soil and solution - Kf - can be used to describe Cu concentrations in the sorbed and solution phases, defined as Eq. (1):

$$K_f = \frac{Cu_{total}}{Cu_{solution}n} (1)$$
<figure><img src="image_10.png" /></figure>

previous Cu deposition and focused on risks arising from future deposition. To do so, we first reviewed the empirical equations estimating Cu's Kf based on soil properties to highlight generic soil properties governing this partition. From this review, we extracted the best compromise Kf equation to estimate partitioning at the regional scale, which ensures more accurate Kf calculation based on pedo-geochemical data typically recorded in soil surveys, thus mainly available. This allowed us to estimate Cu's Kf values to be used at the European scale based on pedo-geochemical soil surveys without the knowledge of soil Cu total content. We then focused on the current state of the climate and its projected changes over the XXIst century, based on two climate change scenarios. To capture the difficulties in runoff prediction and to disentangle the uncertainties between rainfall prediction and runoff calculations of land surface models, we used a set of simulations provided by the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP). These simulations used different land surface models driven by different climate forcings computed by different climate models. For each scenario and each couple of land surface model and climate forcing we estimated the LP or AP of each grid points by comparison between the local values of Kf and of runoff to the respective calculated European median.

## 2. Materials and methods

### 2.1. Equations to estimate copper Kf

The rigorous definition of Kf is based on the concentration ratio of sorbed vs solution species (here Cu) at the equilibrium. Yet, for practical reasons of measurement and applicability, Kf is conventionally derived from total Cu and not from sorbed Cu (Degryse et al., 2009). A general form of the Cu partition coefficient between soil and solution - Kf - can be used to describe Cu concentrations in the sorbed and solution phases, defined as Eq. (1):

$$K_f = \frac{Cu_{total}}{Cu_{solution}n} (1)$$
<figure><img src="image_10.png" /></figure>

Where n stands for the variation in binding strength with metal loading (Groenenberg et al., 2010). A low Kf reflects a high proportion of Cu in solution for a given total Cu content of the soil. Kf can vary as a function of different soil parameters (Degryse et al., 2009; Elzinga et al., 1999) and can also be estimated using Eq. (2):

$$log_{10}(K_f) = a_0 + \sum_i a_i log_{10}(X_i) \quad (2)$$

with X_i the different soil parameters and a_i the corresponding associated coefficient to the parameter.

Numerous studies in the literature have attempted the determination of the value of Kf using the Eq. (2) based on statistical relationships between soil pedo-geochemical parameters, Cu in solution and total Cu measurements.

The soil pedo-geochemical parameter X_i and its associated coefficient a_i can differ depending on the study and the data set used for the estimation. For the purposes of this study, Kf is estimated at the European level, so the formula chosen strikes the best balance between the accuracy of the relationship and its applicability on a wide scale. Thus, the equation must:

i) Include only parameters that are measured in large soil surveys

ii) Fit a large range of each soil parameter

iii) Focus on in situ long-term contamination and not on laboratory experiments.

On December 2020 we first ran a bibliographic research on WOS looking for "Cu AND availab*AND soil AND TOPIC function". We then completed this research by examining the references cited in the articles found. We collected the available relationships for estimating Kf on the basis of soil pedo-geochemical characteristics and/or total Cu. We selected only relationships that were based on commonly collected soil pedo-geochemical characteristics, such as soil organic matter (OM) or soil organic carbon (OC), dissolved organic carbon (DOC), cationic exchange capacity (CEC), clay percentage and pH that are the most frequently reported values from large scale soil survey.

2.2 Soil data
<figure><img src="image_11.png" /></figure>

Where n stands for the variation in binding strength with metal loading (Groenenberg et al., 2010). A low Kf reflects a high proportion of Cu in solution for a given total Cu content of the soil. Kf can vary as a function of different soil parameters (Degryse et al., 2009; Elzinga et al., 1999) and can also be estimated using Eq. (2):

$$log_{10}(K_f) = a_0 + \sum_i a_i log_{10}(X_i) \quad (2)$$

with X_i the different soil parameters and a_i the corresponding associated coefficient to the parameter.

Numerous studies in the literature have attempted the determination of the value of Kf using the Eq. (2) based on statistical relationships between soil pedo-geochemical parameters, Cu in solution and total Cu measurements.

The soil pedo-geochemical parameter X_i and its associated coefficient a_i can differ depending on the study and the data set used for the estimation. For the purposes of this study, Kf is estimated at the European level, so the formula chosen strikes the best balance between the accuracy of the relationship and its applicability on a wide scale. Thus, the equation must:

i) Include only parameters that are measured in large soil surveys

ii) Fit a large range of each soil parameter

iii) Focus on in situ long-term contamination and not on laboratory experiments.

On December 2020 we first ran a bibliographic research on WOS looking for "Cu AND availab*AND soil AND TOPIC function". We then completed this research by examining the references cited in the articles found. We collected the available relationships for estimating Kf on the basis of soil pedo-geochemical characteristics and/or total Cu. We selected only relationships that were based on commonly collected soil pedo-geochemical characteristics, such as soil organic matter (OM) or soil organic carbon (OC), dissolved organic carbon (DOC), cationic exchange capacity (CEC), clay percentage and pH that are the most frequently reported values from large scale soil survey.

2.2 Soil data
<figure><img src="image_11.png" /></figure>

This study used European data on various soil parameters, in particular pH and organic carbon (OC), obtained from the Joint Research Centre's (JRC) LUCAS topsoil data. The data set is limited only to the territories of European Union. The aforementioned data set provides information on pH (https://esdac.jrc.ec.europa.eu/content/copper-distribution-topsoils) and OC contents (https://esdac.jrc.ec.europa.eu/content/topsoil-soil-organic-carbon-lucas-eu25). The data has been re-gridded with cdo commands (Schulzweida, 2019) to a spatial resolution of 0.5° (equivalent to approximately 50 km). This was done to match the resolution of the land surface models that were used to estimate the runoff. The resulting runoff data is presented in section 2.3.

## 2.3. Runoff data from land surface models

Runoff is computed in land models from incoming rain- and snow- falls, calculated evapotranspiration, and soil hydrologic capacities. To estimate changes in runoff across century and to reduce uncertainties, we used two typical land-surface schemes models (LSM) – namely ORCHIDEE (Krinner et al., 2005) and LPJmL (Sitch et al., 2003)– and two global circulation models (GCM) providing climate projections – namely IPSL-CM5a (Dufresne et al., 2013) and GDFL-ESM2m (Dunne et al., 2012) – further named CM5a and ESM2m respectively. Our study exploited simulations conducted as part of the Inter-Sectoral Impact Model Intercomparison Project Phase 2b (ISIMIP2b), which supplied simulations of land surface models driven by binding scenarios from 1861 to 2099 (Frieler et al., 2017). Further details of the protocol used can be found at https://www.isimip.org/protocol/#isimip2b. The ISIMIP2b utilizes harmonized climate forcings derived from gridded, daily bias-adjusted climate data of various CMIP5 (5th coupled model intercomparison project) global circulation models (GCMs) (Frieler et al., 2017; Lange, 2016) as well as with the use of global annual atmospheric CO₂ concentration, and harmonized annual land use maps (Goldewijk et al., 2017). The application of bias-corrected climate data ensures that the climate used by the land surface models is consistent with observations over the last 40 years of the historical period. We compared the historical data calculated by the different models with three five-year periods distributed over the XXIst century: the beginning (2001-2005, called historical scenario), a middle scenario (2051-2055) and an end scenario (2091-2095). In order to simulate mid and end century periods (2051-2055 and 2091-2095), we have used two century-scale scenarios called
<figure><img src="image_12.png" /></figure>

This study used European data on various soil parameters, in particular pH and organic carbon (OC), obtained from the Joint Research Centre's (JRC) LUCAS topsoil data. The data set is limited only to the territories of European Union. The aforementioned data set provides information on pH (https://esdac.jrc.ec.europa.eu/content/copper-distribution-topsoils) and OC contents (https://esdac.jrc.ec.europa.eu/content/topsoil-soil-organic-carbon-lucas-eu25). The data has been re-gridded with cdo commands (Schulzweida, 2019) to a spatial resolution of 0.5° (equivalent to approximately 50 km). This was done to match the resolution of the land surface models that were used to estimate the runoff. The resulting runoff data is presented in section 2.3.

## 2.3. Runoff data from land surface models

Runoff is computed in land models from incoming rain- and snow- falls, calculated evapotranspiration, and soil hydrologic capacities. To estimate changes in runoff across century and to reduce uncertainties, we used two typical land-surface schemes models (LSM) – namely ORCHIDEE (Krinner et al., 2005) and LPJmL (Sitch et al., 2003)– and two global circulation models (GCM) providing climate projections – namely IPSL-CM5a (Dufresne et al., 2013) and GDFL-ESM2m (Dunne et al., 2012) – further named CM5a and ESM2m respectively. Our study exploited simulations conducted as part of the Inter-Sectoral Impact Model Intercomparison Project Phase 2b (ISIMIP2b), which supplied simulations of land surface models driven by binding scenarios from 1861 to 2099 (Frieler et al., 2017). Further details of the protocol used can be found at https://www.isimip.org/protocol/#isimip2b. The ISIMIP2b utilizes harmonized climate forcings derived from gridded, daily bias-adjusted climate data of various CMIP5 (5th coupled model intercomparison project) global circulation models (GCMs) (Frieler et al., 2017; Lange, 2016) as well as with the use of global annual atmospheric CO₂ concentration, and harmonized annual land use maps (Goldewijk et al., 2017). The application of bias-corrected climate data ensures that the climate used by the land surface models is consistent with observations over the last 40 years of the historical period. We compared the historical data calculated by the different models with three five-year periods distributed over the XXIst century: the beginning (2001-2005, called historical scenario), a middle scenario (2051-2055) and an end scenario (2091-2095). In order to simulate mid and end century periods (2051-2055 and 2091-2095), we have used two century-scale scenarios called
<figure><img src="image_12.png" /></figure>

Representative Concentration Pathway (RCP). These scenarios have been defined by the Intergovernmental Panel on Climate Change (IPCC) (van Vuuren et al., 2011) and correspond to common socio-economic pathways followed by the world's population. Here, we focused on RCP 2.6, which represents an active reduction of greenhouse gas emissions to comply with the Paris Agreement, and RCP 6.0, which represents more or less business as usual. RCP 2.6 is predicted to produce a radiation forcing of 2.6 W.m⁻², whereas RCP 6.0 would result in a radiation forcing of 6 W.m⁻².

For each combination of LSMs (LPJmL or ORCHIDEE) and GCMs (CM5a or ESM2m), we calculated the mean over 5 years at the beginning (2001 - 2005), mid (2051 - 2055) and end (2091 - 2095) of the XXIst century. The cross scheme of two land surface models and two GCMs enabled us to establish whether estimations of runoff are influenced more by rainfall projection or the representation of soil hydrologic characteristics. When predictions will be driven by soil hydrologic properties, highest differences in runoff predictions are expected between LPmL_CM5a or ESM2m) and ORCHIDEE_CM5a or ESM2m) projections than between LPJmL_CM5a and ORCHIDEE_CM5a or between LMJmL_ESM2m and ORCHIDEE_ESM2m. Contrarily, when predictions will be driven by rainfall projections, highest differences in runoff predictions are expected between LPmL_CM5a and ORCHIDEE_CM5a or between LPJmL_ESM2m and ORCHIDEE_ESM2m projections than between LPJmL_CM5a and LMJmL_ESM2m or between ORCHIDEE_CM5a and ORCHIDEE_ESM2m.

## 2.4. Statistical tests to assess AP and LP areas

AP or LP areas were assessed by comparing the $K_f$ and runoff values of each grid point with its corresponding spatial median. Median runoff was computed for the whole of Europe for each five-year average period studied. LP areas were characterized by low $K_f$ and high runoff, while AP areas were characterized by the opposite (see Eq. (3a) and (3b)). We employed the classical approach described by Reimann et al., (2005) by classifying as outliers values higher than the median +2 x (median average deviation) (MAD) or lower than the median -2 x MAD. Rather than excluding data points as outliers, we identified data points with unusually high or low values, later referred as anomalies. Thus, we used a threshold lower than 2 MAD for deviation definition and chose to fix a 1 MAD
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

Representative Concentration Pathway (RCP). These scenarios have been defined by the Intergovernmental Panel on Climate Change (IPCC) (van Vuuren et al., 2011) and correspond to common socio-economic pathways followed by the world's population. Here, we focused on RCP 2.6, which represents an active reduction of greenhouse gas emissions to comply with the Paris Agreement, and RCP 6.0, which represents more or less business as usual. RCP 2.6 is predicted to produce a radiation forcing of 2.6 W.m⁻², whereas RCP 6.0 would result in a radiation forcing of 6 W.m⁻².

For each combination of LSMs (LPJmL or ORCHIDEE) and GCMs (CM5a or ESM2m), we calculated the mean over 5 years at the beginning (2001 - 2005), mid (2051 - 2055) and end (2091 - 2095) of the XXIst century. The cross scheme of two land surface models and two GCMs enabled us to establish whether estimations of runoff are influenced more by rainfall projection or the representation of soil hydrologic characteristics. When predictions will be driven by soil hydrologic properties, highest differences in runoff predictions are expected between LPmL_CM5a or ESM2m) and ORCHIDEE_CM5a or ESM2m) projections than between LPJmL_CM5a and ORCHIDEE_CM5a or between LMJmL_ESM2m and ORCHIDEE_ESM2m. Contrarily, when predictions will be driven by rainfall projections, highest differences in runoff predictions are expected between LPmL_CM5a and ORCHIDEE_CM5a or between LPJmL_ESM2m and ORCHIDEE_ESM2m projections than between LPJmL_CM5a and LMJmL_ESM2m or between ORCHIDEE_CM5a and ORCHIDEE_ESM2m.

## 2.4. Statistical tests to assess AP and LP areas

AP or LP areas were assessed by comparing the $K_f$ and runoff values of each grid point with its corresponding spatial median. Median runoff was computed for the whole of Europe for each five-year average period studied. LP areas were characterized by low $K_f$ and high runoff, while AP areas were characterized by the opposite (see Eq. (3a) and (3b)). We employed the classical approach described by Reimann et al., (2005) by classifying as outliers values higher than the median +2 x (median average deviation) (MAD) or lower than the median -2 x MAD. Rather than excluding data points as outliers, we identified data points with unusually high or low values, later referred as anomalies. Thus, we used a threshold lower than 2 MAD for deviation definition and chose to fix a 1 MAD
<figure><img src="image_13.png" /></figure>

187 threshold. MAD was computed as $median(|x_i| - median(x))$, $x$ being successively runoff and $K_f$ for the $i$ grid points where $K_f$ can be estimated (see Eq. (3a) and (3b)).

189 For each combination of LSM (ORCHIDEE or LPJmL) x GCM (CM5a or ESM2m) and each time period (t=2001-2005; 2051-2055 or 2019-2095) with the two climate change scenarios (RCP 2.6 or RCP 6.0) applied for the periods 2051-2055 and 2091-2095, we have defined LP and AP areas as follows:

192 ☐ Areas with soils exhibiting high potentiality of Cu leaching (LP areas) under 1 MAD threshold (named LP) for a 5 years mean time period $t$ were defined as areas where grid points $i$ have:

194 $\begin{cases} K_f(i) < Median(European K_f) - 1 MAD(European K_f) \\ Runoff(t, i) > Median(European runoff(t) + 1 MAD(European Runoff(t))) \end{cases}$ (3a)

195 ☐ Areas with soils exhibiting low potentiality of leaching corresponding to soils of high Cu accumulation potentiality (AP areas) under 1 MAD threshold (named AP) for a 5 years mean time period $t$ were defined as areas where grid points $i$ have:

198 $\begin{cases} K_f(i) > Median(European K_f) + 1 MAD(European K_f) \\ Runoff(t, i) < Median(European runoff(t) - 1 MAD(European Runoff(t))) \end{cases}$ (3b)

199
200 The benefits of this approach is that it is not affected by the set of coefficients chosen to compute $K_f$, and it removes the absolute nature of the values, but focus on the highest (and lowest) values. We used R 4.1.2 (R Core Team, 2021) to compute anomalies and perform the figures.

203
204 3. Results

205 3.1. $K_f$ estimations at the European scale

206
<figure><img src="image_15.png" /></figure>

187 threshold. MAD was computed as $median(|x_i| - median(x))$, $x$ being successively runoff and $K_f$ for the $i$ grid points where $K_f$ can be estimated (see Eq. (3a) and (3b)).

189 For each combination of LSM (ORCHIDEE or LPJmL) x GCM (CM5a or ESM2m) and each time period (t=2001-2005; 2051-2055 or 2019-2095) with the two climate change scenarios (RCP 2.6 or RCP 6.0) applied for the periods 2051-2055 and 2091-2095, we have defined LP and AP areas as follows:

192 ☐ Areas with soils exhibiting high potentiality of Cu leaching (LP areas) under 1 MAD threshold (named LP) for a 5 years mean time period $t$ were defined as areas where grid points $i$ have:

194 $$ \begin{cases} K_f(i) < Median (European K_f) - 1 MAD (European K_f) \\ Runoff(t, i) > Median (European runoff (t) + 1 MAD (European Runoff(t)) \end{cases} \quad (3a) $$

195 ☐ Areas with soils exhibiting low potentiality of leaching corresponding to soils of high Cu accumulation potentiality (AP areas) under 1 MAD threshold (named AP) for a 5 years mean time period $t$ were defined as areas where grid points $i$ have:

198 $$ \begin{cases} K_f(i) > Median (European K_f) + 1 MAD (European K_f) \\ Runoff(t, i) < Median (European runoff (t) - 1 MAD (European Runoff(t)) \end{cases} \quad (3b) $$

199
200 The benefits of this approach is that it is not affected by the set of coefficients chosen to compute $K_f$, and it removes the absolute nature of the values, but focus on the highest (and lowest) values. We used R 4.1.2 (R Core Team, 2021) to compute anomalies and perform the figures.

203
204 3. Results
205 3.1. $K_f$ estimations at the European scale

206
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

The empirical equations extracted from our literature review to estimate K<sub>f</sub> are given in Table 1. We collected 15 equations allowing to calculate K<sub>f</sub> as the coefficient of partition between total Cu and Cu in solution. Among these equations, pH was found the more decisive factor in K<sub>f</sub> estimation (8/15 relationships). Indeed, K<sub>f</sub> is positively correlated to pH with a partial slope for pH around 0.3 for four of these eight relationships so that the more alkaline the soil is, the highest the ratio total Cu/Cu in solution is. Soil organic matter (OM) or OC is less often a parameter in the K<sub>f</sub> equations (4/15 relationships) but, when present, partial slope for OM/OC is higher than that for pH. Three of the 4 papers concerned found a positive relationship between OM and K<sub>f</sub> while Mondaca et al., (2015) found a negative partial slope for soil OM or dissolved OC (Table 1, Eq. (12d)). However, this Eq. (12d) concerns chilies soils and includes a positive partial slope for the CEC. The CEC value can be viewed as the sum of clay and soil OM contents, so that the over whole partial slope of OM is compensated in that particular situation.
<figure><img src="image_16.png" /></figure>

The empirical equations extracted from our literature review to estimate K<sub>f</sub> are given in Table 1. We collected 15 equations allowing to calculate K<sub>f</sub> as the coefficient of partition between total Cu and Cu in solution. Among these equations, pH was found the more decisive factor in K<sub>f</sub> estimation (8/15 relationships). Indeed, K<sub>f</sub> is positively correlated to pH with a partial slope for pH around 0.3 for four of these eight relationships so that the more alkaline the soil is, the highest the ratio total Cu/Cu in solution is. Soil organic matter (OM) or OC is less often a parameter in the K<sub>f</sub> equations (4/15 relationships) but, when present, partial slope for OM/OC is higher than that for pH. Three of the 4 papers concerned found a positive relationship between OM and K<sub>f</sub> while Mondaca et al., (2015) found a negative partial slope for soil OM or dissolved OC (Table 1, Eq. (12d)). However, this Eq. (12d) concerns chilies soils and includes a positive partial slope for the CEC. The CEC value can be viewed as the sum of clay and soil OM contents, so that the over whole partial slope of OM is compensated in that particular situation.
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Table 1.: Transfer functions reviewed from literature to estimate partition coefficient of Cu. R.V stand for response variable and Int. for intercept.
Most studies fitted Kᵢ defined as Kᵢ= [Cu]soil/[Cu]solution n-opt in L.kg⁻¹, Cu.soil or Cu tot in mg.kg⁻¹, DOC (dissolved organic carbon) in mg.L⁻¹, OM (soil organic matter) in %, CEC in cmol.kg⁻¹, standard error around fitted coefficient are reported when indicated in the original article</figcaption></figure>

<table><thead><tr><th>Author</th><th>Eq</th><th>R.V</th><th>Int.</th><th>Log (Cu tot)</th><th>pH</th><th>Log (OM)</th><th>Log (DOC)</th><th>other</th><th>n-opt</th><th>R2</th><th>number of data</th><th>Range Cu tot</th><th>Range OM</th><th>Range DOC</th><th>Range pH</th></tr></thead><tbody><tr><td>(Vulkan et al., 2000)</td><td>4</td><td>Log (K<sub>i</sub>)</td><td>1.74</td><td></td><td>0.34</td><td></td><td>-0.58</td><td></td><td></td><td>0.42</td><td>21</td><td>19-8645</td><td>9.8-69.8</td><td>5.5-8</td><td></td></tr><tr><td>(Sauvé et al., 2000)</td><td>5a</td><td>Log (K<sub>i</sub>)</td><td>1.49 ±0.13</td><td></td><td>0.27 ±0.02</td><td></td><td></td><td></td><td></td><td>0.29</td><td>447</td><td>6.8-82850</td><td></td><td></td><td></td></tr><tr><td>(Sauvé et al., 2000)</td><td>5b</td><td>Log (K<sub>i</sub>)</td><td>1.75 ±0.12</td><td></td><td>0.21 ±0.02</td><td>0.51 ±0.06</td><td></td><td></td><td></td><td>0.42</td><td>353</td><td>6.8-82850</td><td></td><td></td><td></td></tr><tr><td>(Degryse et al., 2009)</td><td>6a</td><td>Log (K<sub>i</sub>)</td><td>0.6</td><td></td><td>0.37</td><td></td><td></td><td>0.65 log (CEC %)</td><td></td><td>0.34</td><td>129</td><td></td><td></td><td></td><td></td></tr><tr><td>(Degryse et al., 2009)</td><td>6b</td><td>Log (K<sub>i</sub>)</td><td>0.45</td><td></td><td>0.34</td><td></td><td></td><td></td><td></td><td>0.44</td><td>128</td><td></td><td></td><td></td><td></td></tr><tr><td>(Unamuno et al., 2009)</td><td>7a</td><td>Log (K<sub>i</sub>)</td><td>1.95</td><td></td><td>0.16</td><td></td><td></td><td></td><td></td><td>0.15</td><td>29</td><td>18-10389</td><td></td><td></td><td></td></tr><tr><td>(Unamuno et al., 2009)</td><td>7b</td><td>Log (K<sub>i</sub>)</td><td>2.383</td><td>0.46</td><td></td><td></td><td></td><td></td><td></td><td>0.61</td><td>29</td><td>18-10389</td><td></td><td></td><td></td></tr><tr><td>(Unamuno et al., 2009)</td><td>7c</td><td>Log (K<sub>i</sub>)</td><td>1.99</td><td>0.42</td><td>0.06</td><td></td><td></td><td></td><td></td><td>0.63</td><td>29</td><td>18-10389</td><td></td><td></td><td></td></tr><tr><td>(Groenenberg et al., 2010)</td><td>8a</td><td>Log (K<sub>i</sub>)</td><td>2.26</td><td></td><td>0.89</td><td>0.9</td><td></td><td></td><td></td><td>0.85</td><td>216</td><td>0.1-326</td><td>2-97.8</td><td>3.3-8.3</td><td></td></tr><tr><td>(Ivezić et al., 2012)</td><td>9a</td><td>Log (K<sub>i</sub>)</td><td>3.98</td><td></td><td></td><td>0.48</td><td>-0.59</td><td></td><td></td><td>0.5</td><td>74</td><td>5.7-141</td><td>0.9-10.2</td><td>4.3-8.1</td><td></td></tr><tr><td>(Mondaca et al., 2015)</td><td>10a</td><td>Log (K<sub>i</sub>)</td><td>1.05</td><td>0.7</td><td></td><td>-1.06</td><td></td><td></td><td></td><td>0.46</td><td>86</td><td>56-4441</td><td>12.0-62</td><td>6.2-7.8</td><td></td></tr><tr><td>(Mondaca et al., 2015)</td><td>10b</td><td>Log (K<sub>i</sub>)</td><td>2.88</td><td>0.41</td><td></td><td></td><td>-1.03</td><td></td><td></td><td>0.77</td><td>86</td><td>56-4441</td><td>12.0-62</td><td>6.2-7.8</td><td></td></tr><tr><td>(Li et al., 2017)</td><td>11a</td><td>Log (K<sub>i</sub>)</td><td>3.12</td><td>0.47</td><td></td><td></td><td>-0.66</td><td></td><td></td><td>0.28</td><td>34</td><td></td><td></td><td></td><td></td></tr><tr><td>(Li et al., 2017)</td><td>11b</td><td>Log (K<sub>i</sub>)</td><td>2.179</td><td>-0.45 * log (Cu solution) µmol.L<sup>-1</sup></td><td></td><td></td><td></td><td></td><td></td><td>0.42</td><td>34</td><td></td><td></td><td></td><td></td></tr><tr><td>(Li et al., 2017)</td><td>11c</td><td>Log (K<sub>i</sub>)</td><td>2.59</td><td>0.617</td><td></td><td></td><td>-1.55</td><td></td><td></td><td>0.88</td><td>20</td><td></td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Table 1.: Transfer functions reviewed from literature to estimate partition coefficient of Cu. R.V stand for response variable and Int. for intercept. Most studies fitted Kᵢ defined as Kᵢ= [Cu]soil/[Cu]solution n-opt in L.kg⁻¹, Cu.soil or Cu tot in mg.kg⁻¹, DOC (dissolved organic carbon) in mg.L⁻¹, OM (soil organic matter) in %, CEC in cmol.kg⁻¹, standard error around fitted coefficient are reported when indicated in the original article</figcaption></figure>

<table><thead><tr><th>Author</th><th>Eq</th><th>R.V</th><th>Int.</th><th>Log (Cu tot)</th><th>pH</th><th>Log (OM)</th><th>Log (DOC)</th><th>other</th><th>n-opt</th><th>R2</th><th>number of data</th><th>Range Cu tot</th><th>Range OM</th><th>Range DOC</th><th>Range pH</th></tr></thead><tbody><tr><td>(Vulkan et al., 2000)</td><td>4</td><td>Log (K<sub>i</sub>)</td><td>1.74</td><td></td><td>0.34</td><td></td><td>-0.58</td><td></td><td></td><td>0.42</td><td>21</td><td>19-8645</td><td>9.8-69.8</td><td>5.5-8</td><td></td></tr><tr><td>(Sauvé et al., 2000)</td><td>5a</td><td>Log (K<sub>i</sub>)</td><td>1.49 ±0.13</td><td></td><td>0.27 ±0.02</td><td></td><td></td><td></td><td></td><td>0.29</td><td>447</td><td>6.8-82850</td><td></td><td></td><td></td></tr><tr><td>(Sauvé et al., 2000)</td><td>5b</td><td>Log (K<sub>i</sub>)</td><td>1.75 ±0.12</td><td></td><td>0.21 ±0.02</td><td>0.51 ±0.06</td><td></td><td></td><td></td><td>0.42</td><td>353</td><td>6.8-82850</td><td></td><td></td><td></td></tr><tr><td>(Degryse et al., 2009)</td><td>6a</td><td>Log (K<sub>i</sub>)</td><td>0.6</td><td></td><td>0.37</td><td></td><td></td><td>0.65 log (CEC %)</td><td></td><td>0.34</td><td>129</td><td></td><td></td><td></td><td></td></tr><tr><td>(Degryse et al., 2009)</td><td>6b</td><td>Log (K<sub>i</sub>)</td><td>0.45</td><td></td><td>0.34</td><td></td><td></td><td></td><td></td><td>0.44</td><td>128</td><td></td><td></td><td></td><td></td></tr><tr><td>(Unamuno et al., 2009)</td><td>7a</td><td>Log (K<sub>i</sub>)</td><td>1.95</td><td></td><td>0.16</td><td></td><td></td><td></td><td></td><td>0.15</td><td>29</td><td>18-10389</td><td></td><td></td><td></td></tr><tr><td>(Unamuno et al., 2009)</td><td>7b</td><td>Log (K<sub>i</sub>)</td><td>2.383</td><td>0.46</td><td></td><td></td><td></td><td></td><td></td><td>0.61</td><td>29</td><td>18-10389</td><td></td><td></td><td></td></tr><tr><td>(Unamuno et al., 2009)</td><td>7c</td><td>Log (K<sub>i</sub>)</td><td>1.99</td><td>0.42</td><td>0.06</td><td></td><td></td><td></td><td></td><td>0.63</td><td>29</td><td>18-10389</td><td></td><td></td><td></td></tr><tr><td>(Groenenberg et al., 2010)</td><td>8a</td><td>Log (K<sub>i</sub>)</td><td>2.26</td><td></td><td>0.89</td><td>0.9</td><td></td><td></td><td></td><td>0.85</td><td>216</td><td>0.1-326</td><td>2-97.8</td><td>3.3-8.3</td><td></td></tr><tr><td>(Ivezić et al., 2012)</td><td>9a</td><td>Log (K<sub>i</sub>)</td><td>3.98</td><td></td><td></td><td>0.48</td><td>-0.59</td><td></td><td></td><td>0.5</td><td>74</td><td>5.7-141</td><td>0.9-10.2</td><td>4.3-8.1</td><td></td></tr><tr><td>(Mondaca et al., 2015)</td><td>10a</td><td>Log (K<sub>i</sub>)</td><td>1.05</td><td>0.7</td><td></td><td>-1.06</td><td></td><td></td><td></td><td>0.46</td><td>86</td><td>56-4441</td><td>12.0-62</td><td>6.2-7.8</td><td></td></tr><tr><td>(Mondaca et al., 2015)</td><td>10b</td><td>Log (K<sub>i</sub>)</td><td>2.88</td><td>0.41</td><td></td><td></td><td>-1.03</td><td></td><td></td><td>0.77</td><td>86</td><td>56-4441</td><td>12.0-62</td><td>6.2-7.8</td><td></td></tr><tr><td>(Li et al., 2017)</td><td>11a</td><td>Log (K<sub>i</sub>)</td><td>3.12</td><td>0.47</td><td></td><td></td><td>-0.66</td><td></td><td></td><td>0.28</td><td>34</td><td></td><td></td><td></td><td></td></tr><tr><td>(Li et al., 2017)</td><td>11b</td><td>Log (K<sub>i</sub>)</td><td>2.179</td><td>-0.45 * log (Cu solution) µmol.L<sup>-1</sup></td><td></td><td></td><td></td><td></td><td></td><td>0.42</td><td>34</td><td></td><td></td><td></td><td></td></tr><tr><td>(Li et al., 2017)</td><td>11c</td><td>Log (K<sub>i</sub>)</td><td>2.59</td><td>0.617</td><td></td><td></td><td>-1.55</td><td></td><td></td><td>0.88</td><td>20</td><td></td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_18.png" /></figure>

1

2

3

4

5

6

7

8

9

10

Over the 15 equations, the estimation of Kf according to Sauvé et al., (2000a) with Eq. (5a) or (5b) (Table 1) is the most robust as determined over a wide range of soils (more than 400 points). The estimations are based on a large gradient of in situ total soil Cu concentrations, even though the highest total soil Cu concentration is higher than what was observed in Europe with the JRC's soil survey. The authors proposed two equations based on a compilation of about 400 data points from long-term contaminated samples. One of the equations considers OM values, whereas the other does not due to a lack of information in the gathered data. Finally, due to the well-known importance in OM for binding with Cu, the Eq. (5b) was selected for our application at the Europe scale and Kf was calculated as following:

$$ \log_{10}(K_f) = 1.75 + 0.21 \times pH + 0.51 \log_{10}(OM) $$

11

with Kf in L.Kg⁻¹ and OM being the soil organic matter content calculated as OM = 2 x OC from JRC following

12

Pribyl (2010).

13

Kf values display a range of 4600 to 21500 L.kg⁻¹ with a median value of 9829 L.kg⁻¹. Kf values below 8000 L.kg⁻¹ and above 12000 L.kg⁻¹ respectively represent low and high anomalies for Kf. On the European scale, a heterogeneous distribution can be seen when using equation (5b), as shown in (Fig. 1).

14

15
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

1

2 Over the 15 equations, the estimation of Kf according to Sauvé et al., (2000a) with Eq. (5a) or (5b) (Table 1) is the

3 most robust as determined over a wide range of soils (more than 400 points). The estimations are based on a

4 large gradient of in situ total soil Cu concentrations, even though the highest total soil Cu concentration is higher

5 than what was observed in Europe with the JRC's soil survey. The authors proposed two equations based on a

6 compilation of about 400 data points from long-term contaminated samples. One of the equations considers OM

7 values, whereas the other does not due to a lack of information in the gathered data. Finally, due to the well-

8 known importance in OM for binding with Cu, the Eq. (5b) was selected for our application at the Europe scale

9 and Kf was calculated as following:

$$10 \quad \log_{10}(K_f) = 1.75 + 0.21 \times pH + 0.51 \log_{10}(OM)$$

11 with Kf in L.Kg⁻¹ and OM being the soil organic matter content calculated as OM = 2 x OC from JRC following

12 Pribyl (2010).

13 Kf values display a range of 4600 to 21500 L.kg⁻¹ with a median value of 9829 L.kg⁻¹. Kf values below 8000 L.kg⁻¹

14 and above 12000 L.kg⁻¹ respectively represent low and high anomalies for Kf. On the European scale, a

15 heterogeneous distribution can be seen when using equation (5b), as shown in (Fig. 1).
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 1: Map of log₁₀(Kf) in Europe at 0.5° following Eq. (5b) applied to soil Cu contents. White pixels correspond to pixel without OC measurement, then no Kf calcul.</figcaption></figure>

Beyond the EU's administrative borders (e.g. Switzerland and Norway), in certain mountain areas there is a lack of OC data which isn't supplied by the JRC. Cu partitioning in soil solution is low around the Mediterranean, UK, Baltic and Nordic regions with high Kf (>12000 L.kg⁻¹). This accounts for 29.9 % of the grid cells, where deposited Cu can thus accumulate in soils. On the contrary, high partition of Cu into soil solution can be found in 20.1 % of the grid cells where values of Kf are low (<8000 L.kg⁻¹), thus providing soils with a tendency of acting as a source of copper for other ecosystems, depending on the runoff. This occurs for instance near Portugal and Poland.

3.2. Modelling potential Cu leaching and accumulation in European soils at the beginning of the century (2001-2005)
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 1: Map of log₁₀(Kf) in Europe at 0.5° following Eq. (5b) applied to soil Cu contents. White pixels correspond to pixel without OC measurement, then no Kf calcul.</figcaption></figure>

Beyond the EU's administrative borders (e.g. Switzerland and Norway), in certain mountain areas there is a lack of OC data which isn't supplied by the JRC. Cu partitioning in soil solution is low around the Mediterranean, UK, Baltic and Nordic regions with high Kf (>12000 L.kg⁻¹). This accounts for 29.9 % of the grid cells, where deposited Cu can thus accumulate in soils. On the contrary, high partition of Cu into soil solution can be found in 20.1 % of the grid cells where values of Kf are low (<8000 L.kg⁻¹), thus providing soils with a tendency of acting as a source of copper for other ecosystems, depending on the runoff. This occurs for instance near Portugal and Poland.

3.2. Modelling potential Cu leaching and accumulation in European soils at the beginning of the century (2001-2005)
<figure><img src="image_22.png" /></figure>

Over the two LSMs x 2 GCMs, the runoff values during the 2001-2005 period varied between 0 (LPJmL_CM5a and LPJmL_ESM2m) and 5.4 mm.day⁻¹ (LPJmL_CM5a). The mean runoff value over the two LSMs x 2 GCMs is 1.1 (± 0.1 standard deviation) mm.day⁻¹ (data shown in Fig S1). For this period, the 1MAD threshold gives rather similar low and high runoff anomalies between couples of LSMs x GCMs, below 0.6, 0.6, 0.7, 0.6 mm.day⁻¹ and above 1.3, 1.2, 1.3 and 1.1 mm.day⁻¹ respectively for ORCHIDEE_CM5a, ORCHIDEE_ESM2m, LPJmL_CM5a and LPJmL_ESM2m. In addition, respectively 21.7, 22.1, 20.2 and 21.1 % of the grid cells are low runoff anomalies and 28.2, 27.9, 29.8 and 28.9 % of the grid cells are high runoff anomalies (see Table S1).

Fig. 2 represents the LP and AP areas for the 2001-2005 period and for the different combinations of LSMs and GCMs. The amount of grid cells concerned by LP and AP areas varied among the LSMs x GCMs combinations (Fig. 3 with the historical scenario and Table S1). However, spatial patterns are well conserved with more similarities between the same LSM than between the same GCM. Globally, LP areas are located mostly in Northern Portugal with scattered points around France, Germany and Scandinavia while AP areas are mostly found in South East of Spain, South-Adriatic coast of Italy and scattered points in Hungary. But, with the ORCHIDEE LSM, AP areas in South Spain are larger, and LP areas in France and East Europe are more scattered than with the LPJmL LSM.
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

Over the two LSMs x 2 GCMs, the runoff values during the 2001-2005 period varied between 0 (LPJmL_CM5a and LPJmL_ESM2m) and 5.4 mm.day⁻¹ (LPJmL_CM5a). The mean runoff value over the two LSMs x 2 GCMs is 1.1 (± 0.1 standard deviation) mm.day⁻¹ (data shown in Fig S1). For this period, the 1MAD threshold gives rather similar low and high runoff anomalies between couples of LSMs x GCMs, below 0.6, 0.6, 0.7, 0.6 mm.day⁻¹ and above 1.3, 1.2, 1.3 and 1.1 mm.day⁻¹ respectively for ORCHIDEE_CM5a, ORCHIDEE_ESM2m, LPJmL_CM5a and LPJmL_ESM2m. In addition, respectively 21.7, 22.1, 20.2 and 21.1 % of the grid cells are low runoff anomalies and 28.2, 27.9, 29.8 and 28.9 % of the grid cells are high runoff anomalies (see Table S1).

Fig. 2 represents the LP and AP areas for the 2001-2005 period and for the different combinations of LSMs and GCMs. The amount of grid cells concerned by LP and AP areas varied among the LSMs x GCMs combinations (Fig. 3 with the historical scenario and Table S1). However, spatial patterns are well conserved with more similarities between the same LSM than between the same GCM. Globally, LP areas are located mostly in Northern Portugal with scattered points around France, Germany and Scandinavia while AP areas are mostly found in South East of Spain, South-Adriatic coast of Italy and scattered points in Hungary. But, with the ORCHIDEE LSM, AP areas in South Spain are larger, and LP areas in France and East Europe are more scattered than with the LPJmL LSM.
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 2: Areas of potential for Cu leaching (LP) and accumulation (AP) over the historical (2001-2005) period for the combinations of land surface scheme (ORCHIDEE in (a), (b); LPJmL in (c), (d)) and climate forcing (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no $K_f$ calcul.</figcaption></figure>

Over the four combinations of LSMs and GCM, 6.4 ± 0.1 % (median, median deviation) of the grid cells are concerned by LP (Fig. 3 (a)) and 6.7 ± 1.1 % by AP (Fig. 3(b)). Areas concerned by LP are almost equal between all LSMs x GCMs even if ESM2m forcing leads to slightly less areas concerned by LP than CM5a. Much more AP areas are predicted by ORCHIDEE LSM. LPJmL_CM5a combination has the smallest percentage of the grid cells concerned by AP with 5.5 %, while ORCHIDEE_CM5A has the largest percentage with 8.0 % (Fig. 3(b)).

3.3. Modelling the evolutions of the LP areas over the century according to the different RCPs
<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 2: Areas of potential for Cu leaching (LP) and accumulation (AP) over the historical (2001-2005) period for the combinations of land surface scheme (ORCHIDEE in (a), (b); LPJmL in (c), (d)) and climate forcing (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no Kf calcul.</figcaption></figure>

Over the four combinations of LSMs and GCM, 6.4 ± 0.1 % (median, median deviation) of the grid cells are concerned by LP (Fig. 3 (a)) and 6.7 ± 1.1 % by AP (Fig. 3(b)). Areas concerned by LP are almost equal between all LSMs x GCMs even if ESM2m forcing leads to slightly less areas concerned by LP than CM5a. Much more AP areas are predicted by ORCHIDEE LSM. LPJmL_CM5a combination has the smallest percentage of the grid cells concerned by AP with 5.5 %, while ORCHIDEE_CM5A has the largest percentage with 8.0 % (Fig. 3(b)).

3.3. Modelling the evolutions of the LP areas over the century according to the different RCPs
<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

For the two chosen climate change scenarios, median runoffs are expected to increase over the century for the 2 LSMs x 2 GCMs combinations. For the mid 2051-2055 period, predicted runoff is 1.1 ± 0.1 mm.day⁻¹ with RCP 2.6 and RCP 6.0 (mean, standard deviation of the 2 LSMs x 2 GCMs over the 5 years), (see Fig. S2 for RCP 2.6 and Fig. S4 for RCP 6.0). For the end 2091-2095 period, predicted runoff is also 1.1 ± 0.1 mm.day⁻¹ with RCP 2.6 but 1.0 ± 0.1 mm.day⁻¹ with RCP 6.0 (mean, standard deviation of the 2 LSMs x 2 GCMs over the 5 years), (Fig. S3 for RCP 2.6 and Fig. S5 for RCP 6.0). Table S1 shows that the amount of grid cells defined as high anomalies for runoff tends to decrease by the end of the century while the amount of grid cells defined as low anomalies for runoff tends to increase. However, tendencies for the 2051-2055 period are variable with in some cases an increase or a decrease in percentage by comparison with the previous or subsequent periods (see Table S1). Furthermore, among the different periods of climate change scenarios, the ratio of LP areas in percentage over areas of high anomalies for runoff is not constant (see Table S1).

The evolution of areas in Europe concerned by LP for the different climate scenarios and the different LSMs x GCMs combinations over the century is presented in percentage in Fig. 3(a). Compared to the historical values and whatever the scenario, the median percentage of grid cells concerned by LP in 2091-2095 decreases by 1.2 ± 0.3 percentage points (median, median deviation) for RCP 2.6 and by 2.1 ± 0.5 percentage points for RCP 6.0. Hence, at the end of the century, percentage of surfaces concerned by LP are 5.3 ± 0.3 % (median, median deviation) for RCP 2.6 and 4.3 ± 0.6 % for RCP 6.0. Estimations of areas concerned by LP are relatively similar for all the time period and climate change scenarios and for all LSMs x GCMs except ORCHIDEE_ESM2m that always predicted smallest percentage of areas concerned by LP. Indeed, for ORCHIDEE_ESM2m the percentage of areas concerned by LP are from 59% (RCP 6.0 2091-2095) to 79 % (RCP 6.0 2051-2055) smallest than the median percentage of surfaces concerned by LP (see Fig. 3(a)).
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

For the two chosen climate change scenarios, median runoffs are expected to increase over the century for the 2 LSMs x 2 GCMs combinations. For the mid 2051-2055 period, predicted runoff is 1.1 ± 0.1 mm.day⁻¹ with RCP 2.6 and RCP 6.0 (mean, standard deviation of the 2 LSMs x 2 GCMs over the 5 years), (see Fig. S2 for RCP 2.6 and Fig. S4 for RCP 6.0). For the end 2091-2095 period, predicted runoff is also 1.1 ± 0.1 mm.day⁻¹ with RCP 2.6 but 1.0 ± 0.1 mm.day⁻¹ with RCP 6.0 (mean, standard deviation of the 2 LSMs x 2 GCMs over the 5 years), (Fig. S3 for RCP 2.6 and Fig. S5 for RCP 6.0). Table S1 shows that the amount of grid cells defined as high anomalies for runoff tends to decrease by the end of the century while the amount of grid cells defined as low anomalies for runoff tends to increase. However, tendencies for the 2051-2055 period are variable with in some cases an increase or a decrease in percentage by comparison with the previous or subsequent periods (see Table S1). Furthermore, among the different periods of climate change scenarios, the ratio of LP areas in percentage over areas of high anomalies for runoff is not constant (see Table S1).

The evolution of areas in Europe concerned by LP for the different climate scenarios and the different LSMs x GCMs combinations over the century is presented in percentage in Fig. 3(a). Compared to the historical values and whatever the scenario, the median percentage of grid cells concerned by LP in 2091-2095 decreases by 1.2 ± 0.3 percentage points (median, median deviation) for RCP 2.6 and by 2.1 ± 0.5 percentage points for RCP 6.0. Hence, at the end of the century, percentage of surfaces concerned by LP are 5.3 ± 0.3 % (median, median deviation) for RCP 2.6 and 4.3 ± 0.6 % for RCP 6.0. Estimations of areas concerned by LP are relatively similar for all the time period and climate change scenarios and for all LSMs x GCMs except ORCHIDEE_ESM2m that always predicted smallest percentage of areas concerned by LP. Indeed, for ORCHIDEE_ESM2m the percentage of areas concerned by LP are from 59% (RCP 6.0 2091-2095) to 79 % (RCP 6.0 2051-2055) smallest than the median percentage of surfaces concerned by LP (see Fig. 3(a)).
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

Fig. 3: Percentage of the grid cells concerned by Cu LP (a) and AP (b) for the different scenarios (historical=2001-2005, RCP 2.6 horizon 2050 and 2090 and RCP 6.0 horizon 2090). The 4 combinations of the 2 LSMs (ORCHIDEE in green and LPJmL in orange) and the 2 climate forcings (CM5a fill bars and ESM2m dashed bar) as well as the median (purple) of the 4 models and the median deviation (bar) are plotted.

The median evolution of LP during the century depends on the climate change scenario. With RCP 2.6, the median percentage of grid cells concerned by LP varied more between the historical scenario and the 2051-2055 one (-0.8 ± 0.4 percentage points, median, median deviation) than between the 2051-2055 and the 2091-2095 periods (-0.4 ± 0.3 percentage points). On the contrary, with RCP 6.0, the median percentage of grid cells concerned by LP areas decreases less from the historical scenario to the 2051-2055 one (-0.3 ± 0.2 percentage points, median, median deviation), than between the 2051-2055 and 2091-2095 periods (-2.0 ± 0.2 percentage points), see Fig. 3 (a). Furthermore, with RCP 2.6, estimations give 5.5 ± 0.5 % of the grid cells concerned by LP in 2051-2055 and 6.2 ± 0.2% with RCP 6.0, which is similar to the 2001-2005 estimate.

For all LSMs and GCMs and the two RCPs, LP areas mostly concern Portugal, north Germany and Scandinavia. In terms of LP risks, the combinations of GCMs and climate change scenarios mostly affect the quantity of dispersed
<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 3: Percentage of the grid cells concerned by Cu LP (a) and AP (b) for the different scenarios (historical=2001-2005, RCP 2.6 horizon 2050 and 2090 and RCP 6.0 horizon 2090). The 4 combinations of the 2 LSMs (ORCHIDEE in green and LPJmL in orange) and the 2 climate forcings (CM5a fill bars and ESM2m dashed bar) as well as the median (purple) of the 4 models and the median deviation (bar) are plotted.</figcaption></figure>

The median evolution of LP during the century depends on the climate change scenario. With RCP 2.6, the median percentage of grid cells concerned by LP varied more between the historical scenario and the 2051-2055 one (-0.8 ± 0.4 percentage points, median, median deviation) than between the 2051-2055 and the 2091-2095 periods (-0.4 ± 0.3 percentage points). On the contrary, with RCP 6.0, the median percentage of grid cells concerned by LP areas decreases less from the historical scenario to the 2051-2055 one (-0.3 ± 0.2 percentage points, median, median deviation), than between the 2051-2055 and 2091-2095 periods (-2.0 ± 0.2 percentage points), see Fig. 3 (a). Furthermore, with RCP 2.6, estimations give 5.5 ± 0.5 % of the grid cells concerned by LP in 2051-2055 and 6.2 ± 0.2% with RCP 6.0, which is similar to the 2001-2005 estimate.

For all LSMs and GCMs and the two RCPs, LP areas mostly concern Portugal, north Germany and Scandinavia. In terms of LP risks, the combinations of GCMs and climate change scenarios mostly affect the quantity of dispersed
<figure><img src="image_29.png" /></figure>

spots in East Europe and the southern extend of Portugal. By 2050, the decrease in LP areas mostly concerns center of France, south of Portugal and north of Germany (Fig. 4 for the RCP 2.6 and Fig. 6 for the RCP 6.0). By 2090, the decrease in LP areas mostly concerns south of Portugal (Fig. 5 for the RCP 2.6 and Fig. 7 for the RCP 6.0).

<figure><img src="image_30.png" /><figcaption>Fig. 4: Areas of potential for Cu leaching (LP) and accumulation (AP) over the RCP2.6 2051-2055 period for the different combinations of land surface schemes (ORCHIDEE in (a), (b) ; LPJmL in (c), (d)) and climate forcings (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no Kf calcul</figcaption></figure>
<figure><img src="image_32.png" /></figure>

spots in East Europe and the southern extend of Portugal. By 2050, the decrease in LP areas mostly concerns center of France, south of Portugal and north of Germany (Fig. 4 for the RCP 2.6 and Fig. 6 for the RCP 6.0). By 2090, the decrease in LP areas mostly concerns south of Portugal (Fig. 5 for the RCP 2.6 and Fig. 7 for the RCP 6.0).

<figure><img src="image_33.png" /></figure>

Fig. 4: Areas of potential for Cu leaching (LP) and accumulation (AP) over the RCP2.6 2051-2055 period for the different combinations of land surface schemes (ORCHIDEE in (a), (b) ; LPJmL in (c), (d)) and climate forcings (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no Kf calcul
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /><figcaption>Fig. 5: Areas of potential for Cu leaching (LP) and accumulation (AP) over the RCP 2.6 2091-2095 period for the different combinations of land surface schemes (ORCHIDEE in (a), (b); LPJmL in (c), (d)) and climate forcings (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no $K_f$ calcul.</figcaption></figure>

### 3.4. Modelling the evolutions of the AP areas over the century according to the different RCPs

The evolution of areas in Europe concerned by AP for the different climate scenarios and the different LSMs x GCMs over the century is presented in percentage in Fig. 3(b). To the end of the century (2091-2095) and for the two climate change scenarios, the percentage of grid cells concerned by AP increases for all LSMs x GCMs except for ORCHIDEE_CM5a with RCP 2.6. AP area increases are highly variable between LSMs x GCMs, with a
<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /><figcaption>Fig. 5: Areas of potential for Cu leaching (LP) and accumulation (AP) over the RCP 2.6 2091-2095 period for the different combinations of land surface schemes (ORCHIDEE in (a), (b) ; LPJmL in (c), (d)) and climate forcings (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no Kf calcul.</figcaption></figure>

3.4. Modelling the evolutions of the AP areas over the century according to the different RCPs

The evolution of areas in Europe concerned by AP for the different climate scenarios and the different LSMs
x GCMs over the century is presented in percentage in Fig. 3(b). To the end of the century (2091-2095) and for
the two climate change scenarios, the percentage of grid cells concerned by AP increases for all LSMs x GCMs
except for ORCHIDEE_CM5a with RCP 2.6. AP area increases are highly variable between LSMs x GCMs, with a
<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

108 smaller increase between historical period and 2091-2095 for RCP 2.6 than for RCP 6.0.

109 With RCP 2.6, and for all LSMs x GCMs, the percentage of surfaces concerned by AP increases between the
110 historical scenario and the mid-one (2051-2055). Between the mid and the 2091-2095 scenarios, the percentage
111 of grid cells concerned by AP increases for LSMs_ESM2m and decreases for LSMs_CM5a (see Fig. 3 (b)).

112 With RCP 6.0, the percentage of areas concerned by AP increases for all LSM x GCM except with ORCHIDEE_CM5a
113 between the historical scenario and the mid-one, and for all LSM x GCM combinations between the mid- and the
114 2091-2095 scenarios.

115 For all LSMs X GCMs and the two RCPs, AP areas are found in Sicilia, East Europe and South Spain. However, the
density and extent of the AP areas in these regions varied between LSMs x GCMs and climate change scenarios
116 (Fig. 4 and 5 for the RCP 2.6 respectively by 2050 and by 2090 and Fig. 6 and 7 for the RCP 6.0 respectively by
117 2050 and by 2090). Over the century, we found new AP areas in East Europe and Greece.
118
<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

108 smaller increase between historical period and 2091-2095 for RCP 2.6 than for RCP 6.0.

109 With RCP 2.6, and for all LSMs x GCMs, the percentage of surfaces concerned by AP increases between the
110 historical scenario and the mid-one (2051-2055). Between the mid and the 2091-2095 scenarios, the percentage
111 of grid cells concerned by AP increases for LSMs_ESM2m and decreases for LSMs_CM5a (see Fig. 3 (b)).

112 With RCP 6.0, the percentage of areas concerned by AP increases for all LSM x GCM except with ORCHIDEE_CM5a
113 between the historical scenario and the mid-one, and for all LSM x GCM combinations between the mid- and the
114 2091-2095 scenarios.

115 For all LSMs X GCMs and the two RCPs, AP areas are found in Sicilia, East Europe and South Spain. However, the
density and extent of the AP areas in these regions varied between LSMs x GCMs and climate change scenarios
116 (Fig. 4 and 5 for the RCP 2.6 respectively by 2050 and by 2090 and Fig. 6 and 7 for the RCP 6.0 respectively by
117 2050 and by 2090). Over the century, we found new AP areas in East Europe and Greece.
118
<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /><figcaption>Fig. 6: Area of potential of leaching (LP) and accumulation (AP) over the RCP 6.0 2051-2055 period for the different combination of land surface scheme (ORCHIDEE in (a), (b) ; LPJmL in (c), (d)) and climate forcings (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no Kf calcul.</figcaption></figure>

119

120

121

122

123
<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /><figcaption>Fig. 6: Area of potential of leaching (LP) and accumulation (AP) over the RCP 6.0 2051-2055 period for the different combination of land surface scheme (ORCHIDEE in (a), (b) ; LPJmL in (c), (d)) and climate forcings (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no $K_f$ calcul.</figcaption></figure>
<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /><figcaption>Fig. 7: Areas of Cu potential for leaching (LP) and accumulation (AP) potential over the RCP 6.0 2091-2095 period for the different combinations of land surface schemes (ORCHIDEE in (a), (b) ; LPJmL in (c), (d)) and climate forcings (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no Kf calcul.</figcaption></figure>

Finally, over all LSMs x GCMs and climate change scenarios, the extent of areas of LP and AP in each region rather depends on GCM than on LSM, with more similarities between ORCHIDEE_GCM (sub figures (a) and (b) in Fig. 2, 4, 5, 6, 7) and LPJmL_GCM (sub figures (c) and (d) in Figs. 2, 4, 5, 6, 7) than between LSM_CM5a (sub figures (a) and (c) in Figs. 2, 4, 5, 6, 7) and LSM_ESM2m (sub figures (b) and (d) in Figs. 2, 4, 5, 6, 7).

124

125

126

127

128

129

130

131

132

133
<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /><figcaption>Fig. 7: Areas of Cu potential for leaching (LP) and accumulation (AP) potential over the RCP 6.0 2091-2095 period for the different combinations of land surface schemes (ORCHIDEE in (a), (b) ; LPJmL in (c), (d)) and climate forcings (CM5a in (a), (c) and ESM2m in (b), (d)). White pixels correspond to pixel without OC measurement, then no Kf calcul.</figcaption></figure>

Finally, over all LSMs x GCMs and climate change scenarios, the extent of areas of LP and AP in each region rather depends on GCM than on LSM, with more similarities between ORCHIDEE_GCM (sub figures (a) and (b) in Fig. 2, 4, 5, 6, 7) and LPJmL_GCM (sub figures (c) and (d) in Figs. 2, 4, 5, 6, 7) than between LSM_CM5a (sub figures (a) and (c) in Figs. 2, 4, 5, 6, 7) and LSM_ESM2m (sub figures (b) and (d) in Figs. 2, 4, 5, 6, 7).

124

125

126

127

128

129

130

131

132

133
<figure><img src="image_40.png" /></figure>

134

135 4. Discussion

136

137 4.1. Modelling the role of copper sink or source with time for contaminated soils

138 This study aims at identifying potential leaching soil areas for Cu over Europe in order to identify locations where
139 soil may play a role in the Cu transfer from soil to aquatic ecosystems. To estimate the proportion of Cu reaching
140 soil solution, we chose to focus on the partitioning coefficient that considers soil properties. This specific choice
141 of Kf coefficient rather than considering only the soil total Cu contents was made because Cu in solution is not
142 strictly correlated with total Cu, nor with other single soil properties as for instance pH and soil OM which are
143 both known to affect Cu partitioning and mobility. Thus, taking into account the variability of soil properties at
144 the European scale, the spatial distribution of Cu in solution was shown to be different from the spatial
145 distribution of total Cu (Sereni et al., 2022). Moreover, data on Cu in solution at large scales are not available
146 making impossible the estimation of AP or LP areas without using the Kr. Finally, the use of partition coefficient
147 allowed us to estimate risk areas without considering total soil Cu temporal variability and with the hypothesis
148 that pedological soil characteristics will not change at the time scale studied. This is a strong implicit assumption
149 but needed at that stage. Indeed even though some soil OM projections are available (Varney et al., 2022) to our
150 knowledge, future projections of pH values at European scale due to climate change are not available limiting our
151 capacities to calculate a time-dependent Kr. Furthermore, together with rainfall and soil moisture changes, climate
152 change is expected to also induce higher temperatures and shorter winters, so that a shift in cultures toward
153 North is expected (Hannah et al., 2013). Therefore, areas with currently low total soil Cu levels may potentially
154 experience a rise in Cu inputs from fungicides, which may subsequently be transported through freshwater
155 systems. Thus, the estimations of LP and AP as emphasized here, can be used to identify high-risk regions and
156 anticipate total content modifications that could occur with an eventual change in anthropogenic activities.
157 Indeed, land management changes due to land use changes or regulation changes may affect the use of Cu in
158 agriculture in the future with potential consequences on Cu leaching.
<figure><img src="image_43.png" /></figure>

134

135 4. Discussion

136

137 4.1. Modelling the role of copper sink or source with time for contaminated soils

138 This study aims at identifying potential leaching soil areas for Cu over Europe in order to identify locations where
139 soil may play a role in the Cu transfer from soil to aquatic ecosystems. To estimate the proportion of Cu reaching
140 soil solution, we chose to focus on the partitioning coefficient that considers soil properties. This specific choice
141 of Kf coefficient rather than considering only the soil total Cu contents was made because Cu in solution is not
142 strictly correlated with total Cu, nor with other single soil properties as for instance pH and soil OM which are
143 both known to affect Cu partitioning and mobility. Thus, taking into account the variability of soil properties at
144 the European scale, the spatial distribution of Cu in solution was shown to be different from the spatial
145 distribution of total Cu (Sereni et al., 2022). Moreover, data on Cu in solution at large scales are not available
146 making impossible the estimation of AP or LP areas without using the Kr. Finally, the use of partition coefficient
147 allowed us to estimate risk areas without considering total soil Cu temporal variability and with the hypothesis
148 that pedological soil characteristics will not change at the time scale studied. This is a strong implicit assumption
149 but needed at that stage. Indeed even though some soil OM projections are available (Varney et al., 2022) to our
150 knowledge, future projections of pH values at European scale due to climate change are not available limiting our
151 capacities to calculate a time-dependent Kr. Furthermore, together with rainfall and soil moisture changes, climate
152 change is expected to also induce higher temperatures and shorter winters, so that a shift in cultures toward
153 North is expected (Hannah et al., 2013). Therefore, areas with currently low total soil Cu levels may potentially
154 experience a rise in Cu inputs from fungicides, which may subsequently be transported through freshwater
155 systems. Thus, the estimations of LP and AP as emphasized here, can be used to identify high-risk regions and
156 anticipate total content modifications that could occur with an eventual change in anthropogenic activities.
157 Indeed, land management changes due to land use changes or regulation changes may affect the use of Cu in
158 agriculture in the future with potential consequences on Cu leaching.
<figure><img src="image_41.png" /></figure>

As a first step, the study conducted here could be used to highlight areas needing regulations to lower Cu input thresholds. Indeed, the evolutions of the LP (and AP) areas we noticed are not only the reflection of the general runoff evolution or of the current Cu risk but also underline areas of interest when combining risk linked to soil contamination and climate change. For instance, in Eastern Europe low K<sub>f</sub> and high runoff result in Cu LP areas with soils tending to act as source of Cu for the other ecosystems. However, in these cases, low amounts of total soil Cu contents (Ballabio et al., 2018) limit the amount of Cu exports. In parallel, in Italy, we found high AP areas whatever the modelling for at least one studied period and one RCP examined. In these vineyard regions, annual Cu inputs are high, resulting in Cu accumulation in soil surface horizons with soils acting as a sink of Cu contamination. These high total Cu concentrations could further enter the food web (García-Esparza et al., 2006) or be exported with soil particles (Imfeld et al., 2020) due to rain erosion (El Azzi et al., 2013). Highly erosive storm events predicted to increase during the next decades in Europe are another risk factor for freshwater contamination even in AP areas. Hence, to go further on, localization of areas with exogenous risks of Cu dissemination have to be identified to reinforce the predictions, e.g. by coupling studies of leaching potential as the one we conducted here with erosion risk studies (Ballabio et al., 2017) and with retention pond localization areas.

174

175 4.2. Temporal evolution of data and scope of the modelling analysis

176 To reduce intra and inter annual variability the modelling conducted here focused on 5-years means, thus aimed
177 at smoothing seasonal variability of runoff or of Cu inputs. The K<sub>f</sub> we calculated was not a dynamic value since we
178 did not make hypothesis about the temporal evolution of soil organic carbon or pH. Furthermore, K<sub>f</sub> is defined on
179 the assumption that there is equilibrium between the solid and solution phases. This means that the amount of
180 Cu in solution estimated by this method may be less than that present immediately after Cu application and before
181 equilibrium is reached (McBride et al., 1997). Despite all, our results showed a good agreement between the four
182 LSMs x GCMs in their projection of amount of grid cells concerned by both LP and AP, validating the use of their
183 median to perform projections in the absence of in situ validation. The scope of our predictions had limits that
<figure><img src="image_44.png" /></figure>

As a first step, the study conducted here could be used to highlight areas needing regulations to lower Cu input thresholds. Indeed, the evolutions of the LP (and AP) areas we noticed are not only the reflection of the general runoff evolution or of the current Cu risk but also underline areas of interest when combining risk linked to soil contamination and climate change. For instance, in Eastern Europe low K<sub>f</sub> and high runoff result in Cu LP areas with soils tending to act as source of Cu for the other ecosystems. However, in these cases, low amounts of total soil Cu contents (Ballabio et al., 2018) limit the amount of Cu exports. In parallel, in Italy, we found high AP areas whatever the modelling for at least one studied period and one RCP examined. In these vineyard regions, annual Cu inputs are high, resulting in Cu accumulation in soil surface horizons with soils acting as a sink of Cu contamination. These high total Cu concentrations could further enter the food web (García-Esparza et al., 2006) or be exported with soil particles (Imfeld et al., 2020) due to rain erosion (El Azzi et al., 2013). Highly erosive storm events predicted to increase during the next decades in Europe are another risk factor for freshwater contamination even in AP areas. Hence, to go further on, localization of areas with exogenous risks of Cu dissemination have to be identified to reinforce the predictions, e.g. by coupling studies of leaching potential as the one we conducted here with erosion risk studies (Ballabio et al., 2017) and with retention pond localization areas.

174

175 4.2. Temporal evolution of data and scope of the modelling analysis

176 To reduce intra and inter annual variability the modelling conducted here focused on 5-years means, thus aimed
177 at smoothing seasonal variability of runoff or of Cu inputs. The K<sub>f</sub> we calculated was not a dynamic value since we
178 did not make hypothesis about the temporal evolution of soil organic carbon or pH. Furthermore, K<sub>f</sub> is defined on
179 the assumption that there is equilibrium between the solid and solution phases. This means that the amount of
180 Cu in solution estimated by this method may be less than that present immediately after Cu application and before
181 equilibrium is reached (McBride et al., 1997). Despite all, our results showed a good agreement between the four
182 LSMs x GCMs in their projection of amount of grid cells concerned by both LP and AP, validating the use of their
183 median to perform projections in the absence of in situ validation. The scope of our predictions had limits that
<figure><img src="image_42.png" /></figure>

rely on the difficulties to predict whether rain- and snow-falls and runoff will evolve in terms of intensity and frequency, even if alternations of drying and rewetting events may affect Cu partitioning between phases (Christensen and Christensen, 2003; Han et al., 2001). To gain field reality at the scale of territory for example, modelling will require to account for the time periods of year with higher rain- and snowfalls amounts coinciding with periods of Cu use, for instance in agriculture and vineyards (Banas et al., 2010; Ribolzi et al., 2002). Indeed, if intense rainfall occurs close to Cu fungicide applications, a larger Cu amount than expected may be exported through runoff. Thus, regional soil Cu budgets require the use of temporal model, which accounts for the regular inputs and outputs of Cu from vegetation and runoff. Finally, the identification of the areas with high risks of soil Cu leaching or accumulation we made in this study can be viewed as a first step for the risk evolution assessment of Cu contamination useful for land management or Cu-fertilizer applications regulations.

## 5. Conclusion

Our approach to assess European areas with a potential to accumulate or leach copper from soils was not straightforward but included several steps. We focused first on the means to calculate Cu partitioning. By reviewing existing Cu $K_f$'s equations we pointed out pH and soil OM contents as important determinants and more precisely that the OM partial effect was larger than the pH one. Then, using the European maps of soil characteristic data, we computed the map of $K_f$ at the 0.5° scale, highlighting areas with high risk to leach or to accumulate Cu for a given soil's overall content or upcoming. The estimation of LP and AP areas for current and future soil runoffs under two RCPs with couples of two GCMs x 2 LSMs was thereafter performed by comparing anomalies for both $K_f$ and runoffs. Interestingly, our first result showed that the variations in the number of LP and AP grid points was not only due to variations in the runoff intensities distribution but also to their localization. Indeed, the ratio of percentage of areas of LP or AP over areas of high or low anomalies was not constant during the century. At the beginning of the XXIst century our study showed that comparable amounts of grid cells are concerned by LP and by AP (between [6.2% - 6.4%] and between [5.5% - 8.0%] respectively). During the century, AP areas were found to increase for all the LSMs x GCMs and the two RCPs. On the contrary, for the two RCPs and
<figure><img src="image_45.png" /></figure>

rely on the difficulties to predict whether rain- and snow-falls and runoff will evolve in terms of intensity and frequency, even if alternations of drying and rewetting events may affect Cu partitioning between phases (Christensen and Christensen, 2003; Han et al., 2001). To gain field reality at the scale of territory for example, modelling will require to account for the time periods of year with higher rain- and snowfalls amounts coinciding with periods of Cu use, for instance in agriculture and vineyards (Banas et al., 2010; Ribolzi et al., 2002). Indeed, if intense rainfall occurs close to Cu fungicide applications, a larger Cu amount than expected may be exported through runoff. Thus, regional soil Cu budgets require the use of temporal model, which accounts for the regular inputs and outputs of Cu from vegetation and runoff. Finally, the identification of the areas with high risks of soil Cu leaching or accumulation we made in this study can be viewed as a first step for the risk evolution assessment of Cu contamination useful for land management or Cu-fertilizer applications regulations.

194

195 5. Conclusion

196 Our approach to assess European areas with a potential to accumulate or leach copper from soils was not straightforward but included several steps. We focused first on the means to calculate Cu partitioning. By reviewing existing Cu K<sub>f</sub>'s equations we pointed out pH and soil OM contents as important determinants and more precisely that the OM partial effect was larger than the pH one. Then, using the European maps of soil characteristic data, we computed the map of K<sub>f</sub> at the 0.5° scale, highlighting areas with high risk to leach or to accumulate Cu for a given soil's overall content or upcoming. The estimation of LP and AP areas for current and future soil runoffs under two RCPs with couples of two GCMs x 2 LSMs was thereafter performed by comparing anomalies for both K<sub>f</sub> and runoffs. Interestingly, our first result showed that the variations in the number of LP and AP grid points was not only due to variations in the runoff intensities distribution but also to their localization. Indeed, the ratio of percentage of areas of LP or AP over areas of high or low anomalies was not constant during the century. At the beginning of the XXI<sup>st</sup> century our study showed that comparable amounts of grid cells are concerned by LP and by AP (between [6.2% - 6.4%] and between [5.5% - 8.0%] respectively). During the century, AP areas were found to increase for all the LSMs x GCMs and the two RCPs. On the contrary, for the two RCPs and
<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

209 three over the four LSMs x GCMs, LP areas were found to decrease during the century compared to the current estimation. Projections for 2090 with RCP 2.6 considering median (and median deviation), indicated 5.3 ± 0.3% of the grid cells concerned by LP areas and 7.9 ± 1.3% by AP areas. Projections for 2090 with RCP 6.0 showed a slightly smallest amount of grid cells concerned by LP and a highest by AP with respectively 4.3 ± 0.6% and 10.3 ± 0.7% of the grid cells (median, median deviation). Surprisingly, the total amount of grid cells concerned by the two risks of AP and LP is rather similar between the two climate change scenarios with estimation between 13.2 ± 1.3 and 14.6 ± 1.3%. This was due, however, to opposite trends in the evolution of LP and AP areas. Their relative proportions and period of main variations differed with most of the evolution between historical period and 2051-2055 for RCP 2.6 and between 2051-2055 and 2091-2095 for RCP 6.0. Finally, we showed that even if the number of grid points identified with LP and AP may varied between LSMs x GCMs models and climate change scenarios, their localizations are roughly conserved, emphasizing the necessity to precise monitoring in Cu application on these areas.

221

222 Code availability:

223 The code can be provided upon request.

224

225 Data availability:

226 The data can be provided upon request

227

228 Credit authorships contribution statement:

229 Laura Sereni: Methodology, Formal analysis, Data processing, Writing original draft.

230 Julie-Mai Paris: Formal analysis, Initial data processing, Writing original draft.

231 Isabelle Lamy: Methodology Conceptualization, Writing review and editing, Supervision, Funding acquisition
<figure><img src="image_46.png" /></figure>

209 three over the four LSMs x GCMs, LP areas were found to decrease during the century compared to the current estimation. Projections for 2090 with RCP 2.6 considering median (and median deviation), indicated 5.3 ± 0.3% of the grid cells concerned by LP areas and 7.9 ± 1.3% by AP areas. Projections for 2090 with RCP 6.0 showed a slightly smallest amount of grid cells concerned by LP and a highest by AP with respectively 4.3 ± 0.6% and 10.3 ± 0.7% of the grid cells (median, median deviation). Surprisingly, the total amount of grid cells concerned by the two risks of AP and LP is rather similar between the two climate change scenarios with estimation between 13.2 ± 1.3 and 14.6 ± 1.3%. This was due, however, to opposite trends in the evolution of LP and AP areas. Their relative proportions and period of main variations differed with most of the evolution between historical period and 2051-2055 for RCP 2.6 and between 2051-2055 and 2091-2095 for RCP 6.0. Finally, we showed that even if the number of grid points identified with LP and AP may varied between LSMs x GCMs models and climate change scenarios, their localizations are roughly conserved, emphasizing the necessity to precise monitoring in Cu application on these areas.

221

222 Code availability:

223 The code can be provided upon request.

224

225 Data availability:

226 The data can be provided upon request

227

228 Credit authorships contribution statement:

229 Laura Sereni: Methodology, Formal analysis, Data processing, Writing original draft.

230 Julie-Mai Paris: Formal analysis, Initial data processing, Writing original draft.

231 Isabelle Lamy: Methodology Conceptualization, Writing review and editing, Supervision, Funding acquisition
<figure><img src="image_45.png" /></figure>

232 Bertrand Guenet: Methodology, conceptualization, Writing review and editing, supervision, Project administration,

233

234 Declaration of competing interests

235 The authors declare that they have no known competing financial interests or personal relationships that could

236 have appeared to influence the work reported in this paper.

237

238 Acknowledgments

239 Parts of this study were financially supported by the Labex BASC through the Connexion project. LS thanks the Ecole

240 Normale Supérieure (ENS) for funding her PhD. The authors thank Nathalie de Noblet-Ducoudré for valuable
discussions on this paper.

242

243 Bibliography:

244 El Azzi, D., Viers, J., Guiresse, M., Probst, A., Aubert, D., Caparros, J., Charles, F., Guizien, K. and Probst, J. L.: Origin
and fate of copper in a small Mediterranean vineyard catchment: New insights from combined chemical
extraction and δ65Cu isotopic composition, Sci. Total Environ., 463–464, 91–101,
doi:10.1016/j.scitotenv.2013.05.058, 2013.

248 Babcsányi, I., Chabaux, F., Granet, M., Meite, F., Payraveau, S., Duplay, J. and Imfeld, G.: Copper in soil fractions
and runoff in a vineyard catchment: Insights from copper stable isotopes, Sci. Total Environ., 557–558, 154–162,
doi:10.1016/j.scitotenv.2016.03.037, 2016.

251 Ballabio, C., Borrelli, P., Spinoni, J., Meusburger, K., Michaelides, S., Beguería, S., Klik, A., Petan, S., Janeček, M.,
<figure><img src="image_47.png" /></figure>

232 Bertrand Guenet: Methodology, conceptualization, Writing review and editing, supervision, Project administration,

233

234 Declaration of competing interests

235 The authors declare that they have no known competing financial interests or personal relationships that could

236 have appeared to influence the work reported in this paper.

237

238 Acknowledgments

239 Parts of this study were financially supported by the Labex BASC through the Connexion project. LS thanks the Ecole

240 Normale Supérieure (ENS) for funding her PhD. The authors thank Nathalie de Noblet-Ducoudré for valuable

241 discussions on this paper.

242

243 Bibliography:

244 El Azzi, D., Viers, J., Guiresse, M., Probst, A., Aubert, D., Caparros, J., Charles, F., Guizien, K. and Probst, J. L.: Origin

245 and fate of copper in a small Mediterranean vineyard catchment: New insights from combined chemical

246 extraction and δ65Cu isotopic composition, Sci. Total Environ., 463–464, 91–101,

247 doi:10.1016/j.scitotenv.2013.05.058, 2013.

248 Babcsányi, I., Chabaux, F., Granet, M., Meite, F., Payraudeau, S., Duplay, J. and Imfeld, G.: Copper in soil fractions

249 and runoff in a vineyard catchment: Insights from copper stable isotopes, Sci. Total Environ., 557–558, 154–162,

250 doi:10.1016/j.scitotenv.2016.03.037, 2016.

251 Ballabio, C., Borrelli, P., Spinoni, J., Meusburger, K., Michaelides, S., Beguería, S., Klik, A., Petan, S., Janeček, M.,
<figure><img src="image_46.png" /></figure>

252 Olsen, P., Aalto, J., Lakatos, M., Rymszewicz, A., Dumitrescu, A., Tadić, M. P., Diodato, N., Kostalova, J., Rousseva, S., Banasik, K., Alewell, C. and Panagos, P.: Mapping monthly rainfall erosivity in Europe, Sci. Total Environ., 579, 1298–1315, doi:10.1016/j.scitotenv.2016.11.123, 2017.

255 Ballabio, C., Panagos, P., Lugato, E., Huang, J. H., Orgiazzi, A., Jones, A., Fernández-Ugalde, O., Borrelli, P. and Montanarella, L.: Copper distribution in European topsoils: An assessment based on LUCAS soil survey, Sci. Total Environ., 636, 282–298, doi:10.1016/j.scitotenv.2018.04.268, 2018.

258 Banas, D., Marin, B., Skraber, S., Chopin, E. I. B. and Zanella, A.: Copper mobilization affected by weather conditions in a stormwater detention system receiving runoff waters from vineyard soils (Champagne, France), Environ. Pollut., 158(2), 476–482, doi:10.1016/j.envpol.2009.08.034, 2010.

261 Broos, K., Warne, M. S. J., Heemsbergen, D. A., Stevens, D., Barnes, M. B., Correll, R. L. and McLaughlin, M. J.: Soil factors controlling the toxicity of copper and zinc to microbial processes in Australian soils, Environ. Toxicol. Chem., 26(4), 583–590, doi:10.1897/06-302R.1, 2007.

264 Christensen, J. H. and Christensen, O. B.: Severe summertime flooding in Europe, Nature, 421(6925), 805–806, doi:10.1038/421805a, 2003.

266 Chu, H., Wei, J., Qiu, J., Li, Q. and Wang, G.: Identification of the impact of climate change and human activities on rainfall-runoff relationship variation in the Three-River Headwaters region, Ecol. Indic., 106, doi:10.1016/j.ecolind.2019.105516, 2019.

269 Degryse, F., Smolders, E. and Parker, D. R.: Partitioning of metals (Cd, Co, Cu, Ni, Pb, Zn) in soils: concepts, methodologies, prediction and applications - a review, Eur. J. Soil Sci., 60(4), 590–612, doi:10.1111/j.1365-2389.2009.01142.x, 2009.

272 Douville, H., Raghavan, K., Renwick, J., Allan, R. P., Arias, P. A., Barlow, M., Cerezo-Mota, R., Cherchi, T., Gan, A. Y.,
<figure><img src="image_48.png" /></figure>

252 Olsen, P., Aalto, J., Lakatos, M., Rymszewicz, A., Dumitrescu, A., Tadić, M. P., Diodato, N., Kostalova, J., Rousseva, S., Banasik, K., Alewell, C. and Panagos, P.: Mapping monthly rainfall erosivity in Europe, Sci. Total Environ., 579, 1298–1315, doi:10.1016/j.scitotenv.2016.11.123, 2017.

255 Ballabio, C., Panagos, P., Lugato, E., Huang, J. H., Orgiazzi, A., Jones, A., Fernández-Ugalde, O., Borrelli, P. and Montanarella, L.: Copper distribution in European topsoils: An assessment based on LUCAS soil survey, Sci. Total Environ., 636, 282–298, doi:10.1016/j.scitotenv.2018.04.268, 2018.

258 Banas, D., Marin, B., Skraber, S., Chopin, E. I. B. and Zanella, A.: Copper mobilization affected by weather conditions in a stormwater detention system receiving runoff waters from vineyard soils (Champagne, France), Environ. Pollut., 158(2), 476–482, doi:10.1016/j.envpol.2009.08.034, 2010.

261 Broos, K., Warne, M. S. J., Heemsbergen, D. A., Stevens, D., Barnes, M. B., Correll, R. L. and McLaughlin, M. J.: Soil factors controlling the toxicity of copper and zinc to microbial processes in Australian soils, Environ. Toxicol. Chem., 26(4), 583–590, doi:10.1897/06-302R.1, 2007.

264 Christensen, J. H. and Christensen, O. B.: Severe summertime flooding in Europe, Nature, 421(6925), 805–806, doi:10.1038/421805a, 2003.

266 Chu, H., Wei, J., Qiu, J., Li, Q. and Wang, G.: Identification of the impact of climate change and human activities on rainfall-runoff relationship variation in the Three-River Headwaters region, Ecol. Indic., 106, doi:10.1016/j.ecolind.2019.105516, 2019.

269 Degryse, F., Smolders, E. and Parker, D. R.: Partitioning of metals (Cd, Co, Cu, Ni, Pb, Zn) in soils: concepts, methodologies, prediction and applications - a review, Eur. J. Soil Sci., 60(4), 590–612, doi:10.1111/j.1365-2389.2009.01142.x, 2009.

272 Douville, H., Raghavan, K., Renwick, J., Allan, R. P., Arias, P. A., Barlow, M., Cerezo-Mota, R., Cherchi, T., Gan, A. Y.,
<figure><img src="image_47.png" /></figure>

273 Gergis, J., Jiang, D., Khan, A., Pokam Mba, W., Rosenfeld, D., Tierney, J. and Zolina, O.: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, in Fundamental and Applied Climatology, vol. 2, edited by V. Masson-Delmotte, P. Zhai, A. Pirani, S. L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen, L. Goldfarb, M. I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J. B. R. Matthews, T. K. Maycock, T. Waterfield, O. Yelekçi, R. Yu, and B. Zhou, pp. 13–25, cambride university press., 2021.

279 Dufresne, J. L., Foujols, M. A., Denvil, S., Caubel, A., Marti, O., Aumont, O., Balkanski, Y., Bekki, S., Bellenger, H., Benshila, R., Bony, S., Bopp, L., Braconnot, P., Brockmann, P., Cadule, P., Cheruy, F., Codron, F., Cozic, A., Cugnet, D., de Noblet, N., Duvel, J. P., Ethé, C., Fairhead, L., Fichefet, T., Flavoni, S., Friedlingstein, P., Grandpeix, J. Y., Guez, L., Guilyardi, E., Hauglustaine, D., Hourdin, F., Idelkadi, A., Ghattas, J., Joussaume, S., Kageyama, M., Krinner, G., Labetoulle, S., Lahellec, A., Lefebvre, M. P., Lefevre, F., Levy, C., Li, Z. X., Lloyd, J., Lott, F., Madec, G., Mancip, M., Marchand, M., Masson, S., Meurdesoif, Y., Mignot, J., Musat, I., Parouty, S., Polcher, J., Rio, C., Schulz, M., Swingedouw, D., Szopa, S., Talandier, C., Terray, P., Viovy, N. and Vuichard, N.: Climate change projections using the IPSL-CM5 Earth System Model: From CMIP3 to CMIP5, Clim. Dyn., doi:10.1007/s00382-012-1636-1, 2013.

287 Dunne, J. P., John, J. G., Adcroft, A. J., Griffies, S. M., Hallberg, R. W., Shevliakova, E., Stouffer, R. J., Cooke, W., Dunne, K. A., Harrison, M. J., Krasting, J. P., Malyshev, S. L., Milly, P. C. D., Phillipps, P. J., Sentman, L. T., Samuels, B. L., Spelman, M. J., Winton, M., Wittenberg, A. T. and Zadeh, N.: GFDL’s ESM2 global coupled climate-carbon earth system models. Part I: Physical formulation and baseline simulation characteristics, J. Clim., doi:10.1175/JCLI-D-11-00560.1, 2012.

292 Elzinga, E. J., Van Grinsven, J. J. M. and Swartjes, F. A.: General purpose Freundlich isotherms for cadmium, copper and zinc in soils, Eur. J. Soil Sci., 50(1), 139–149, doi:10.1046/j.1365-2389.1999.00220.x, 1999.

294 Flemming, C. A. and Trevors, J. T.: Copper toxicity and chemistry in the environment: a review, Water. Air. Soil Pollut., 44(1–2), 143–158, doi:10.1007/BF00228784, 1989.
<figure><img src="image_49.png" /></figure>

273 Gergis, J., Jiang, D., Khan, A., Pokam Mba, W., Rosenfeld, D., Tierney, J. and Zolina, O.: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, in Fundamental and Applied Climatology, vol. 2, edited by V. Masson-Delmotte, P. Zhai, A. Pirani, S. L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen, L. Goldfarb, M. I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J. B. R. Matthews, T. K. Maycock, T. Waterfield, O. Yelekçi, R. Yu, and B. Zhou, pp. 13–25, cambride university press., 2021.

279 Dufresne, J. L., Foujols, M. A., Denvil, S., Caubel, A., Marti, O., Aumont, O., Balkanski, Y., Bekki, S., Bellenger, H., Benshila, R., Bony, S., Bopp, L., Braconnot, P., Brockmann, P., Cadule, P., Cheruy, F., Codron, F., Cozic, A., Cugnet, D., de Noblet, N., Duvel, J. P., Ethé, C., Fairhead, L., Fichefet, T., Flavoni, S., Friedlingstein, P., Grandpeix, J. Y., Guez, L., Guilyardi, E., Hauglustaine, D., Hourdin, F., Idelkadi, A., Ghattas, J., Joussaume, S., Kageyama, M., Krinner, G., Labetoulle, S., Lahellec, A., Lefebvre, M. P., Lefevre, F., Levy, C., Li, Z. X., Lloyd, J., Lott, F., Madec, G., Mancip, M., Marchand, M., Masson, S., Meurdesoif, Y., Mignot, J., Musat, I., Parouty, S., Polcher, J., Rio, C., Schulz, M., Swingedouw, D., Szopa, S., Talandier, C., Terray, P., Viovy, N. and Vuichard, N.: Climate change projections using the IPSL-CM5 Earth System Model: From CMIP3 to CMIP5, Clim. Dyn., doi:10.1007/s00382-012-1636-1, 2013.

287 Dunne, J. P., John, J. G., Adcroft, A. J., Griffies, S. M., Hallberg, R. W., Shevliakova, E., Stouffer, R. J., Cooke, W., Dunne, K. A., Harrison, M. J., Krasting, J. P., Malyshev, S. L., Milly, P. C. D., Phillipps, P. J., Sentman, L. T., Samuels, B. L., Spelman, M. J., Winton, M., Wittenberg, A. T. and Zadeh, N.: GFDL’s ESM2 global coupled climate-carbon earth system models. Part I: Physical formulation and baseline simulation characteristics, J. Clim., doi:10.1175/JCLI-D-11-00560.1, 2012.

292 Elzinga, E. J., Van Grinsven, J. J. M. and Swartjes, F. A.: General purpose Freundlich isotherms for cadmium, copper and zinc in soils, Eur. J. Soil Sci., 50(1), 139–149, doi:10.1046/j.1365-2389.1999.00220.x, 1999.

294 Flemming, C. A. and Trevors, J. T.: Copper toxicity and chemistry in the environment: a review, Water. Air. Soil Pollut., 44(1–2), 143–158, doi:10.1007/BF00228784, 1989.
<figure><img src="image_48.png" /></figure>

296 Frieler, K., Lange, S., Piontek, F., Reyer, C. P. O., Schewe, J., Warszawski, L., Zhao, F., Chini, L., Denvil, S., Emanuel, K., Geiger, T., Halladay, K., Hurtt, G., Mengel, M., Murakami, D., Ostberg, S., Popp, A., Riva, R., Stevanovic, M., SuzGBRi, T., Volkholz, J., Burke, E., Ciais, P., Ebi, K., Eddy, T. D., Elliott, J., Galbraith, E., Gosling, S. N., Hattermann, F., Hickler, T., Hinkel, J., Hof, C., Huber, V., Jägermeyr, J., Krysanova, V., Marcé, R., Müller Schmied, H., Mouratiadou, I., Pierson, D., Tittensor, D. P., Vautard, R., Van Vliet, M., Biber, M. F., Betts, R. A., Leon Bodirsky, B., Deryng, D., Frolking, S., Jones, C. D., Lotze, H. K., Lotze-Campen, H., Sahajpal, R., Thonicke, K., Tian, H. and Yamagata, Y.: Assessing the impacts of 1.5°C global warming - Simulation protocol of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP2b), Geosci. Model Dev., 10(12), 4321–4345, doi:10.5194/gmd-10-4321-2017, 2017.

305 García-Esparza, M. A., Capri, E., Pirzadeh, P. and Trevisan, M.: Copper content of grape and wine from Italian farms, Food Addit. Contam., 23(3), 274–280, doi:10.1080/02652030500429117, 2006.

307 Giller, K. E., Witter, E. and Mcgrath, S. P.: Toxicity of heavy metals to microorganisms and microbial processes in agricultural soils: A review, Soil Biol. Biochem., 30(10–11), 1389–1414, doi:10.1016/S0038-0717(97)00270-8, 1998.

310 Goldewijk, K. K., Beusen, A., Doelman, J. and Stehfest, E.: Anthropogenic land use estimates for the Holocene - HYDE 3.2, Earth Syst. Sci. Data, 9(2), 927–953, doi:10.5194/essd-9-927-2017, 2017.

312 Groenenberg, J. E., Römkens, P. F. A. M., Comans, R. N. J., Luster, J., Pampura, T., Shotbolt, L., Tipping, E. and De Vries, W.: Transfer functions for solid-solution partitioning of cadmium, copper, nickel, lead and zinc in soils: Derivation of relationships for free metal ion activities and validation with independent data, Eur. J. Soil Sci., 61(1), 58–73, doi:10.1111/j.1365-2389.2009.01201.x, 2010.

316 Han, F. X., Banin, A. and Triplett, G. B.: Redistribution of heavy metals in arid-zone soils under a wetting-drying cycle soil moisture regime, Soil Sci., 166(1), 18–28, doi:10.1097/00010694-200101000-00005, 2001.
<figure><img src="image_50.png" /></figure>

296 Frieler, K., Lange, S., Piontek, F., Reyer, C. P. O., Schewe, J., Warszawski, L., Zhao, F., Chini, L., Denvil, S., Emanuel, K., Geiger, T., Halladay, K., Hurtt, G., Mengel, M., Murakami, D., Ostberg, S., Popp, A., Riva, R., Stevanovic, M., SuzGBRi, T., Volkholz, J., Burke, E., Ciais, P., Ebi, K., Eddy, T. D., Elliott, J., Galbraith, E., Gosling, S. N., Hattermann, F., Hickler, T., Hinkel, J., Hof, C., Huber, V., Jägermeyr, J., Krysanova, V., Marcé, R., Müller Schmied, H., Mouratiadou, I., Pierson, D., Tittensor, D. P., Vautard, R., Van Vliet, M., Biber, M. F., Betts, R. A., Leon Bodirsky, B., Deryng, D., Frolking, S., Jones, C. D., Lotze, H. K., Lotze-Campen, H., Sahajpal, R., Thonicke, K., Tian, H. and Yamagata, Y.: Assessing the impacts of 1.5°C global warming - Simulation protocol of the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP2b), Geosci. Model Dev., 10(12), 4321–4345, doi:10.5194/gmd-10-4321-2017, 2017.

305 García-Esparza, M. A., Capri, E., Pirzadeh, P. and Trevisan, M.: Copper content of grape and wine from Italian farms, Food Addit. Contam., 23(3), 274–280, doi:10.1080/02652030500429117, 2006.

307 Giller, K. E., Witter, E. and Mcgrath, S. P.: Toxicity of heavy metals to microorganisms and microbial processes in agricultural soils: A review, Soil Biol. Biochem., 30(10–11), 1389–1414, doi:10.1016/S0038-0717(97)00270-8, 1998.

310 Goldewijk, K. K., Beusen, A., Doelman, J. and Stehfest, E.: Anthropogenic land use estimates for the Holocene - HYDE 3.2, Earth Syst. Sci. Data, 9(2), 927–953, doi:10.5194/essd-9-927-2017, 2017.

312 Groenenberg, J. E., Römkens, P. F. A. M., Comans, R. N. J., Luster, J., Pampura, T., Shotbolt, L., Tipping, E. and De Vries, W.: Transfer functions for solid-solution partitioning of cadmium, copper, nickel, lead and zinc in soils: Derivation of relationships for free metal ion activities and validation with independent data, Eur. J. Soil Sci., 61(1), 58–73, doi:10.1111/j.1365-2389.2009.01201.x, 2010.

316 Han, F. X., Banin, A. and Triplett, G. B.: Redistribution of heavy metals in arid-zone soils under a wetting-drying cycle soil moisture regime, Soil Sci., 166(1), 18–28, doi:10.1097/00010694-200101000-00005, 2001.
<figure><img src="image_49.png" /></figure>

318 Hannah, L., Roehrdanz, P. R., Ikegami, M., Shepard, A. V., Shaw, M. R., Tabor, G., Zhi, L., Marquet, P. A. and Hijmans, R. J.: Climate change, wine, and conservation, Proc. Natl. Acad. Sci. U. S. A., 110(17), 6907–6912, doi:10.1073/pnas.1210127110, 2013.

321 Imfeld, G., Meite, F., Wiegert, C., Guyot, B., Masbou, J. and Payraveau, S.: Do rainfall characteristics affect the export of copper, zinc and synthetic pesticides in surface runoff from headwater catchments?, Sci. Total Environ., 741, 140437, doi:10.1016/j.scitotenv.2020.140437, 2020.

324 Komárek, M., Čadková, E., Chrastný, V., Bordas, F. and Bollinger, J. C.: Contamination of vineyard soils with fungicides: A review of environmental and toxicological aspects, Environ. Int., 36(1), 138–151, doi:10.1016/j.envint.2009.10.005, 2010.

327 Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S. and Prentice, I. C.: A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system, Global Biogeochem. Cycles, doi:10.1029/2003GB002199, 2005.

330 Lange, S.: Earth2Observe, WFDEI and ERA-Interim data Merged and Bias-corrected for ISIMIP (EWEMBI), GFZ Data Serv., 2016.

332 McBride, M., Sauvé, S. and Hendershot, W.: Solubility control of Cu, Zn, Cd and Pb in contaminated soils, Eur. J. Soil Sci., 48(2), 337–346, doi:10.1111/j.1365-2389.1997.tb00554.x, 1997.

334 Mimikou, M. A., Baltas, E., Varanou, E. and Pantazis, K.: Regional impacts of climate change on water resources quantity and quality indicators, J. Hydrol., 234(1–2), 95–109, doi:10.1016/S0022-1694(00)00244-4, 2000.

336 Mondaca, P., Neaman, A., Sauvé, S., Salgado, E. and Bravo, M.: Solubility, partitioning, and activity of copper-contaminated soils in a semiarid region, J. Plant Nutr. Soil Sci., 178(3), 452–459, doi:10.1002/jpln.201400349, 2015.
<figure><img src="image_51.png" /></figure>

318 Hannah, L., Roehrdanz, P. R., Ikegami, M., Shepard, A. V., Shaw, M. R., Tabor, G., Zhi, L., Marquet, P. A. and Hijmans, R. J.: Climate change, wine, and conservation, Proc. Natl. Acad. Sci. U. S. A., 110(17), 6907–6912, doi:10.1073/pnas.1210127110, 2013.

321 Imfeld, G., Meite, F., Wiegert, C., Guyot, B., Masbou, J. and Payraveau, S.: Do rainfall characteristics affect the export of copper, zinc and synthetic pesticides in surface runoff from headwater catchments?, Sci. Total Environ., 741, 140437, doi:10.1016/j.scitotenv.2020.140437, 2020.

324 Komárek, M., Čadková, E., Chrastný, V., Bordas, F. and Bollinger, J. C.: Contamination of vineyard soils with fungicides: A review of environmental and toxicological aspects, Environ. Int., 36(1), 138–151, doi:10.1016/j.envint.2009.10.005, 2010.

327 Krinner, G., Viovy, N., de Noblet-Ducoudré, N., Ogée, J., Polcher, J., Friedlingstein, P., Ciais, P., Sitch, S. and Prentice, I. C.: A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system, Global Biogeochem. Cycles, doi:10.1029/2003GB002199, 2005.

330 Lange, S.: Earth2Observe, WFDEI and ERA-Interim data Merged and Bias-corrected for ISIMIP (EWEMBI), GFZ Data Serv., 2016.

332 McBride, M., Sauvé, S. and Hendershot, W.: Solubility control of Cu, Zn, Cd and Pb in contaminated soils, Eur. J. Soil Sci., 48(2), 337–346, doi:10.1111/j.1365-2389.1997.tb00554.x, 1997.

334 Mimikou, M. A., Baltas, E., Varanou, E. and Pantazis, K.: Regional impacts of climate change on water resources quantity and quality indicators, J. Hydrol., 234(1–2), 95–109, doi:10.1016/S0022-1694(00)00244-4, 2000.

336 Mondaca, P., Neaman, A., Sauvé, S., Salgado, E. and Bravo, M.: Solubility, partitioning, and activity of copper-contaminated soils in a semiarid region, J. Plant Nutr. Soil Sci., 178(3), 452–459, doi:10.1002/jpln.201400349, 2015.
<figure><img src="image_50.png" /></figure>

339 Noll, M. R.: Trace Elements in Terrestrial Environments., 2003.

340 Pribyl, D. W.: A critical review of the conventional SOC to SOM conversion factor, Geoderma, 156(3–4), 75–83,

341 doi:10.1016/j.geoderma.2010.02.003, 2010.

342 R Core Team: R core team (2021), R A Lang. Environ. Stat. Comput. R Found. Stat. Comput. Vienna, Austria. URL

343 http://www. R-project. org, ISBN 3-900051-07-0, URL http://www.R-project.org/. [online] Available from:

344 http://www.mendeley.com/research/r-language-environment-statistical-computing-

345 96/%5Cnpapers2://publication/uuid/A1207DAB-22D3-4A04-82FB-D4DD5AD57C28, 2021.

346 Reimann, C., Filzmoser, P. and Garrett, R. G.: Background and threshold: Critical comparison of methods of

347 determination, Sci. Total Environ., 346(1–3), 1–16, doi:10.1016/j.scitotenv.2004.11.023, 2005.

348 Ribolzi, O., Valles, V., Gomez, L. and Voltz, M.: Speciation and origin of particulate copper in runoff water from a

349 Mediterranean vineyard catchment, Environ. Pollut., 117(2), 261–271, doi:10.1016/S0269-7491(01)00274-3,

350 2002.

351 Rooney, C. P., Zhao, F. J. and McGrath, S. P.: Soil factors controlling the expression of copper toxicity to plants in a

352 wide range of European soils, Environ. Toxicol. Chem., 25(3), 726–732, doi:10.1897/04-602R.1, 2006.

353 Salminen, R. and Gregorauskiene, V.: Considerations regarding the definition of a geochemical baseline of

354 elements in the surficial materials in areas differing in basic geology, Appl. Geochemistry, doi:10.1016/S0883-

355 2927(99)00077-3, 2000.

356 Sauvé, S., Hendershot, W. and Allen, H. E.: Solid-solution partitioning of metals in contaminated soils: Dependence

357 on pH, total metal burden, and organic matter, Environ. Sci. Technol., 34(7), 1125–1131, doi:10.1021/es9907764,

358 2000.

359 Schulzweida, U.: CDO User Guide, doi:10.5281/zenodo.2558193, 2019.
<figure><img src="image_52.png" /></figure>

339 Noll, M. R.: Trace Elements in Terrestrial Environments., 2003.

340 Pribyl, D. W.: A critical review of the conventional SOC to SOM conversion factor, Geoderma, 156(3–4), 75–83,

341 doi:10.1016/j.geoderma.2010.02.003, 2010.

342 R Core Team: R core team (2021), R A Lang. Environ. Stat. Comput. R Found. Stat. Comput. Vienna, Austria. URL

343 http://www. R-project. org, ISBN 3-900051-07-0, URL http://www.R-project.org/. [online] Available from:

344 http://www.mendeley.com/research/r-language-environment-statistical-computing-

345 96/%5Cnpapers2://publication/uuid/A1207DAB-22D3-4A04-82FB-D4DD5AD57C28, 2021.

346 Reimann, C., Filzmoser, P. and Garrett, R. G.: Background and threshold: Critical comparison of methods of

347 determination, Sci. Total Environ., 346(1–3), 1–16, doi:10.1016/j.scitotenv.2004.11.023, 2005.

348 Ribolzi, O., Valles, V., Gomez, L. and Voltz, M.: Speciation and origin of particulate copper in runoff water from a

349 Mediterranean vineyard catchment, Environ. Pollut., 117(2), 261–271, doi:10.1016/S0269-7491(01)00274-3,

350 2002.

351 Rooney, C. P., Zhao, F. J. and McGrath, S. P.: Soil factors controlling the expression of copper toxicity to plants in a

352 wide range of European soils, Environ. Toxicol. Chem., 25(3), 726–732, doi:10.1897/04-602R.1, 2006.

353 Salminen, R. and Gregorauskiene, V.: Considerations regarding the definition of a geochemical baseline of

354 elements in the surficial materials in areas differing in basic geology, Appl. Geochemistry, doi:10.1016/S0883-

355 2927(99)00077-3, 2000.

356 Sauvé, S., Hendershot, W. and Allen, H. E.: Solid-solution partitioning of metals in contaminated soils: Dependence

357 on pH, total metal burden, and organic matter, Environ. Sci. Technol., 34(7), 1125–1131, doi:10.1021/es9907764,

358 2000.

359 Schulzweida, U.: CDO User Guide, doi:10.5281/zenodo.2558193, 2019.
<figure><img src="image_51.png" /></figure>

360 Sereni, L., Guenet, B. and Lamy, I.: Mapping risks associated with soil copper contamination using availability and bio-availability proxies at the European scale, Environ. Sci. Pollut. Res., (123456789), doi:10.1007/s11356-022-23046-0, 2022.

363 Sitch, S., Smith, B., Prentice, I. C., Arneth, A., Bondeau, A., Cramer, W., Kaplan, J. O., Levis, S., Lucht, W., Sykes, M. T., Thonicke, K. and Venevsky, S.: Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model, Glob. Chang. Biol., doi:10.1046/j.1365-2486.2003.00569.x, 2003.

367 Varney, R. M., Chadburn, S. E., Burke, E. J. and Cox, P. M.: Evaluation of soil carbon simulation in CMIP6 Earth system models, Biogeosciences, doi:10.5194/bg-19-4671-2022, 2022.

369 Vidal, M., Santos, M. J., Abrão, T., Rodríguez, J. and Rigol, A.: Modeling competitive metal sorption in a mineral soil, Geoderma, 149(3–4), 189–198, doi:10.1016/j.geoderma.2008.11.040, 2009.

371 van Vuuren, D. P., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., Hurtt, G. C., Kram, T., Krey, V., Lamarque, J. F., Masui, T., Meinshausen, M., Nakicenovic, N., Smith, S. J. and Rose, S. K.: The representative concentration pathways: An overview, Clim. Change, 109(1), 5–31, doi:10.1007/s10584-011-0148-z, 2011.

374 West, T. S. and Coombs, T. L.: Soil as the Source of Trace Elements [and Discussion], Philos. Trans. R. Soc. Lond. B. Biol. Sci., 294(1071), 19–39 [online] Available from: http://www.jstor.org/stable/2395553, 1981.

376
<figure><img src="image_53.png" /></figure>

360 Sereni, L., Guenet, B. and Lamy, I.: Mapping risks associated with soil copper contamination using availability and bio-availability proxies at the European scale, Environ. Sci. Pollut. Res., (123456789), doi:10.1007/s11356-022-23046-0, 2022.

363 Sitch, S., Smith, B., Prentice, I. C., Arneth, A., Bondeau, A., Cramer, W., Kaplan, J. O., Levis, S., Lucht, W., Sykes, M. T., Thonicke, K. and Venevsky, S.: Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model, Glob. Chang. Biol., doi:10.1046/j.1365-2486.2003.00569.x, 2003.

367 Varney, R. M., Chadburn, S. E., Burke, E. J. and Cox, P. M.: Evaluation of soil carbon simulation in CMIP6 Earth system models, Biogeosciences, doi:10.5194/bg-19-4671-2022, 2022.

369 Vidal, M., Santos, M. J., Abrão, T., Rodríguez, J. and Rigol, A.: Modeling competitive metal sorption in a mineral soil, Geoderma, 149(3–4), 189–198, doi:10.1016/j.geoderma.2008.11.040, 2009.

371 van Vuuren, D. P., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., Hurtt, G. C., Kram, T., Krey, V., Lamarque, J. F., Masui, T., Meinshausen, M., Nakicenovic, N., Smith, S. J. and Rose, S. K.: The representative concentration pathways: An overview, Clim. Change, 109(1), 5–31, doi:10.1007/s10584-011-0148-z, 2011.

374 West, T. S. and Coombs, T. L.: Soil as the Source of Trace Elements [and Discussion], Philos. Trans. R. Soc. Lond. B. Biol. Sci., 294(1071), 19–39 [online] Available from: http://www.jstor.org/stable/2395553, 1981.

376