2.1 Equations to estimate copper K_f

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

where Cu_total is the total Cu content of soil (in mg kg⁻¹), Cu_solution is the Cu content of the soil solution (in mg L⁻¹), and n stands for the variation in binding strength with metal loading (Groenenberg et al., 2010). A low K_f reflects a high proportion of Cu in the solution for a given total Cu content of the soil. K_f 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):

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

Numerous studies in the literature have attempted the determination of the value of K_f using Eq. (2) based on statistical relationships between soil pedo-geochemical parameters, Cu in the 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, K_f is estimated at the European Union level; thus, 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.

  be fitted to 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 conducted bibliographic research on Web of Science (WOS) looking for the following: ‑Cu ” AND ‑availab* ” AND ‑soil ” AND TOPIC function. We then completed this search by examining the references cited in the articles found. We collected the available relationships for estimating K_f 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 – constituting the most frequently reported values from large-scale soil surveys.

2.2 Soil data

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 member states. The aforementioned data set provides information on pH (Panagos et al., 2022; Ballabio et al., 2016; ESDAC – European Commission, 2024; Panagos et al., 2012) and OC contents (Panagos et al., 2022; de Brogniez et al., 2015; ESDAC – European Commission, 2024; Panagos et al., 2012). The data have 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 are presented in Sect. 2.3.

2.3 Runoff data from land surface models

Runoff is computed in land models based on incoming precipitation, calculated evapotranspiration, and soil hydrologic capacities. To estimate changes in soil runoff during the 21st century and to reduce uncertainties, we used two typical land surface scheme models (LSMs) – namely ORCHIDEE (Krinner et al., 2005) and LPJmL (Sitch et al., 2003)– and two global circulation models (GCMs) 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 the period 1861–2099 (Frieler et al., 2017). Further details of the protocol used can be found in the ISIMIP2b documentation (The Inter-Sectoral Impact Model Intercomparison Project, 2021). The ISIMIP2b utilizes harmonized climate forcings derived from gridded, daily bias-adjusted climate data of various CMIP5 (Coupled Model Intercomparison Project 5) global circulation models (GCMs) (Frieler et al., 2017; Lange, 2016), as well as global annual atmospheric CO₂ concentration and harmonized annual land use maps (Klein 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 5-year periods distributed over the 21st century: 2001–2005, called the historical scenario; 2051–2055; and 2091–2095. In order to simulate the 2051–2055 and 2091–2095 periods, we used two century-scale scenarios called representative concentration pathways (RCPs). 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 in 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 for three evenly spaced periods of the 21st century: 2001–2005, 2051–2055, and 2091–2095. The cross scheme of two land surface models and two GCMs enabled us to establish whether estimations of runoff are influenced more by rainfall projections provided by the GCMs or by the representation of soil hydrologic characteristics provided by the LSMs. When predictions are driven by soil hydrologic properties, the highest differences in runoff predictions are expected between couplings of the same LSM with different GCMs (e.g., LPJmL_CM5a is closer to LPJmL_ESM2m than to ORCHIDEE_CM5a). Contrarily, when predictions are driven by rainfall projections, the highest differences in runoff predictions are expected between couplings of the same GCM with different LSMs (e.g., LPJmL_CM5a is closer to ORCHIDEE_CM5a than to LPJml_ESM2m).

2.4 Assessment of AP and LP areas

AP or LP areas were assessed by comparing the K_f and runoff values of each grid cell with their corresponding spatial median. Median runoff was computed for the whole of Europe for each 5-year period studied per model. LP areas were characterized by low K_f and high runoff, while AP areas were characterized by the opposite (see Eqs. 3 and 4). We identified grid cells with unusually high or low values, later referred to as anomalies, as grid cells above or below 1 mean absolute deviation (MAD). MAD was computed as median(|x_i|−median(x)), with x being successive runoff, and K_f was computed for the grid cells where it can be estimated (see Eqs. 3 and 4).

For each combination of LSM (ORCHIDEE or LPJmL) and GCM (CM5a or ESM2m) and for 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:

• Areas with soils exhibiting high potentiality of Cu leaching (LP areas, referred to as LP) under the 1-MAD threshold for a 5-year mean time period t were defined as areas where grid cells i have

  $$\begin{array}{}\text{(3a)}& \begin{array}{rl}{K}_{\mathrm{f}}\left(i\right)& <\mathrm{median}\left(\mathrm{European}{K}_{\mathrm{f}}\right)\\ & -\mathrm{1}\mathrm{MAD}\left(\mathrm{European}{K}_{\mathrm{f}}\right),\\ \text{runoff}\left(t,i\right)& >\mathrm{median}\left(\mathrm{European}\mathrm{runoff}\left(t\right)\right)\\ & +\mathrm{MAD}\left(\mathrm{European}\mathrm{runoff}\left(t\right)\right).\end{array}\end{array}$$

• Areas with soils exhibiting low potentiality of leaching and high Cu accumulation potentiality (AP areas, referred to as AP) under the 1 MAD-threshold for a 5-year mean time period t were defined as areas where grid cells i have

  $$\begin{array}{}\text{(3b)}& \begin{array}{rl}{K}_{\mathrm{f}}\left(i\right)& >\mathrm{median}\left(\mathrm{European}{K}_{\mathrm{f}}\right)\\ & +\mathrm{1}\mathrm{MAD}\left(\mathrm{European}{K}_{\mathrm{f}}\right),\\ \text{runoff}\left(t,i\right)& <\mathrm{median}\left(\mathrm{European}\mathrm{runoff}\left(t\right)\right)\\ & -\mathrm{MAD}\left(\mathrm{European}\mathrm{runoff}\left(t\right)\right).\end{array}\end{array}$$

The benefit of this approach is that anomaly identification is not affected by the set of coefficients chosen to compute K_f, and it removes the absolute nature of the values. However, it focus on the deviation from the median.

We choose to calculate the MAD of each time period to emphasize the spatial variability. Anomaly identification could also be done using the historical runoff as a reference and looking at its change with time. However, when considering the actual rainfall regime as a reference, we consider the fact that the current environmental risk includes the spatial risk variability.

In the next sections, the results of the temporal trends are presented using the median per model and the mean over the four models.

We used R 4.1.2 (R Core Team, 2021) to compute anomalies and to create the figures.

3.1 K_f estimations at the European scale

The empirical equations extracted from our literature review to estimate K_f are given in Table 1. We collected 15 equations allowing us to calculate K_f as the coefficient of partitioning between total Cu and Cu in the solution. Among these equations, pH was found to be the more decisive factor in K_f estimation (8 out of 15 relationships). Indeed, K_f is positively correlated to pH so that the more alkaline the soil is, the higher the ratio of total Cu to Cu in the solution is. Soil organic matter (OM) and OC are parameters in the K_f equations less often (4 out of 15 relationships), but, when present, the partial slope for OM or OC is higher than that for pH, which means that a small variation in soil OM content affects more Cu partitioning than a small variation in pH. Three of the four papers concerned found a positive relationship between OM and K_f, while Mondaca et al. (2015) found a negative partial slope for soil OM or dissolved OC (Table 1, Eq. 10b). However, Eq. (12d) was fitted on arid soils from Chile and includes a positive partial slope for the CEC. The CEC value can be viewed as a proxy for the sum of clay and soil OM contents so that the whole partial slope of OM is compensated for in that particular situation.

Table 1Transfer functions reviewed from the literature to estimate partitioning coefficient of Cu. RV stands for response variable, and Int. stands for intercept. Most studies fitted K_f defined as K_f= [Cu]_soil $/$ [Cu]${}_{\mathrm{solution}}^{n\text{-opt}}$ (in L kg⁻¹), Cu_soil or Cu tot (in mg kg⁻¹), DOC (dissolved organic carbon) (in mg L⁻¹), OM (soil organic matter) (in %), and CEC in (cmol kg⁻¹); the standard errors around fitted coefficients are reported when indicated in the original article.

Download Print Version | Download XLSX

Over the 15 equations, the estimation of K_f according to Sauvé et al. (2000) with Eqs. (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 (Ballabio et al., 2018). Sauvé et al. (2000) 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 of OM for binding with Cu, Eq. (5b) was selected for our application at the Europe scale, and K_f was calculated as follows:

with K_f in liters per kilogram (L kg⁻¹) and with OM being the soil organic matter content calculated as OM = 2 × OC, in accordance with the JRC, following Pribyl (2010).

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

Figure 1Map of log₁₀(K_f) in Europe at 0.5° following Eq. (5b) as applied to soil Cu contents. White pixels correspond to pixels without OC measurements and, consequently, no K_f estimations.

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

3.2 Modeling potential Cu leaching and accumulation in European soils for the historical period (2001–2005)

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

Figure 2 represents the LP and AP areas for the 2001–2005 period and for the different combinations of LSMs and GCMs. The number of grid cells with LP and AP areas varied among the LSM and GCM 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 the southeast of Spain and the southern Adriatic coast of Italy and at scattered points in Hungary. However, with the ORCHIDEE LSM, AP areas in southern Spain are larger, and LP areas in France and eastern Europe are more scattered than with the LPJmL LSM.

Figure 2Areas of potential for Cu leaching (LP) and accumulation (AP) over the historical (2001–2005) period for the combinations of land surface schemes (ORCHIDEE in a, b; LPJmL in c, d) and climate forcings (CM5a in a, c; ESM2m in b, d). White pixels correspond to pixels without OC measurement and, consequently, no K_f estimations.

Figure 3Percentage of the grid cells with 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 four combinations of the two LSMs (ORCHIDEE in green and LPJmL in orange) and the two climate forcings (CM5a fill bars and ESM2m dashed bar), as well as the median (purple) of the four models and the median deviations (bar), are plotted.

Download

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

Figure 4Areas of potential for Cu leaching (LP) and accumulation (AP) over the RCP 2.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; ESM2m in b, d). White pixels correspond to pixels without OC measurements and, consequently, no K_f estimations.

Figure 5Areas 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 pixels without OC measurements and, consequently, no K_f estimations.

3.3 Modeling the changes in the LP areas over the century according to the different RCPs

For the two chosen climate change scenarios, median runoffs per model are expected to increase over the century for the combinations of the two LSMs and the two GCMs. For the 2051–2055 period, predicted runoff is 1.1 ± 0.1 mm d⁻¹ with RCP 2.6 and RCP 6.0 (mean, standard deviation of the two LSMs and the two GCMs over the 5 years) (see Fig. S2 for RCP 2.6 and Fig. S4 for RCP 6.0). For the 2091–2095 period, predicted runoff is also 1.1 ± 0.1 mm d⁻¹ with RCP 2.6 but 1.0 ± 0.1 mm d⁻¹ with RCP 6.0 (mean, standard deviation of the two LSMs and the two GCMs over the 5 years) (Fig. S3 for RCP 2.6 and Fig. S5 for RCP 6.0). Table S1 shows that the number of grid cells defined as high anomalies for runoff tends to decrease by the 2091–2095 period, while the number 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 in comparison to the previous or subsequent periods (see Table S1). Furthermore, among the different periods of climate change scenarios, the ratio of LP areas (in terms of percentage) over areas of high anomalies for runoff is not constant (see Table S1).

The change in the areas in Europe with LP for the different climate scenarios and the different LSM and GCM combinations over the century is presented (in terms of percentage) in Fig. 3a. Compared to the historical values, and whatever the scenario, the median percentage of grid cells with 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, for the 2091–2095 period, the percentage of surfaces with LP is 5.3 ± 0.3 % (median, median deviation) for RCP 2.6 and 4.3 ± 0.6 % for RCP 6.0. Areas where LP was detected are relatively similar for all the time periods and climate change scenarios and for all LSM and GCM combinations, except for ORCHIDEE_ESM2m, which always predicted the smallest percentage of areas with LP. Indeed, for ORCHIDEE_ESM2m, the percentage of areas with LP is from 59 % (RCP 6.0, 2091–2095) to 79 % (RCP 6.0, 2051–2055) smaller than the median percentage of surfaces with LP (see Fig. 3a).

The change in the LP's median during the century depends on the climate change scenario. With RCP 2.6, the median percentage of grid cells with LP varied more between the historical scenario and the 2051–2055 scenario (−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 with LP areas decreases less from the historical scenario to the 2051–2055 scenario (−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. 3a. Furthermore, with RCP 2.6, estimations indicate that 5.5 ± 0.5 % of the grid cells demonstrate LP in 2051–2055, while this increases to 6.2 ± 0.2 % with RCP 6.0, which is similar to the 2001–2005 estimate.

For all LSM and GCM combinations and the two RCPs, LP areas are mostly detected in Portugal, northern Germany, and Scandinavia. In terms of LP risks, the combinations of GCMs and climate change scenarios mostly affect the quantity of dispersed spots in eastern Europe and in the southern regions of Portugal. By 2050, the decrease in LP areas is mostly located in the center of France, south of Portugal, and north of Germany (Fig. 4 for RCP 2.6 and Fig. 6 for RCP 6.0). By 2090, the decrease in LP areas is mostly located in the south of Portugal (Fig. 5 for RCP 2.6 and Fig. 7 for RCP 6.0).

3.4 Modeling the changes in the AP areas over the century according to the different RCPs

The change in the AP areas in Europe with the different climate scenarios and the different LSM and GCM combinations over the century is presented (in terms of percentage) in Fig. 3b. For the 2091–2095 period and for the two climate change scenarios, the percentage of grid cells in which AP is detected increases for all LSM and GCM combinations, except for ORCHIDEE_CM5a with RCP 2.6. AP area increases are highly variable between combinations of LSMs and GCMs, with a smaller increase between the historical period and 2091–2095 for RCP 2.6 than for RCP 6.0.

With RCP 2.6 and for all combinations of LSMs and GCMs, the percentage of grid cells where AP is detected increases between the historical scenario and the 2051–2055 period. Between 2051–2055 and 2091–2095, the percentage of grid cells where AP is detected increases for LSMs_ESM2m and decreases for LSMs_CM5a (see Fig. 3b).

With RCP 6.0, the percentage of areas where AP is detected increases for all LSM × GCM except with ORCHIDEE_CM5a between the historical period and the 2051–2055 period, and for all LSM × GCM combinations between the 2051–2055 and the 2091–2095 period.

For all LSM and GCM combinations and the two RCPs, AP areas are found in Sicily, eastern Europe, and southern Spain. However, the density and extent of the AP areas in these regions varied between LSM and GCM combinations and climate change scenarios (Figs. 4 and 5 for RCP 2.6 for the 2051–2055 and 2091–2095 periods, respectively, and Figs. 6 and 7 for RCP 6.0 for the 20510–2055 and 2091–2095 periods, respectively). Over the century, we found new AP areas in eastern Europe and Greece.

Figure 6Area of potential of leaching (LP) and accumulation (AP) over the RCP 6.0 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; ESM2m in b, d). White pixels correspond to pixels without OC measurements and, consequently, no K_f estimations.

Figure 7Areas of Cu potential for leaching (LP) and accumulation (AP) 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; ESM2m in b, d). White pixels correspond to pixels without OC measurements and, consequently, no K_f estimations.

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

4.1 Modeling soil copper release or storage with time for contaminated soils

This study aims to identify potential leaching areas for Cu over Europe in order to identify locations where soil may play a role in the Cu transfer from soil to aquatic ecosystems. To estimate the proportion of Cu reaching the soil solution, we chose to focus on the partitioning coefficient, which was calculated based on soil properties (pH and OM here) other than total soil Cu. This specific choice of the K_f coefficient rather than considering only the soil total Cu contents was made because Cu in the solution is not strictly correlated with total Cu or with other single soil properties such as pH and soil OM, which are both known to affect Cu partitioning and mobility. Thus, taking into account the variability of soil properties at the European scale, the spatial distribution of Cu in the solution was shown to be different from the spatial distribution of total Cu (Sereni et al., 2022). However, data on Cu in the solution at large scales are not available, making the direct estimation of transport within the soil solution and of AP or LP areas without using the K_f impossible. Finally, the use of the partitioning coefficient allowed us to estimate risk areas without considering total soil Cu temporal variability and with the hypothesis that pedological soil characteristics will not change at the timescale studied. This is a strong implicit assumption, but it is needed at this stage. Indeed, even though some soil OM projections are available (Varney et al., 2022), to our knowledge, future projections of pH values at the European scale due to climate change are not available, limiting our capacity to calculate a time-dependent K_f. In particular, there are large uncertainties about the C stocks that may change as a result of climate change and dedicated policies for increasing the C stocks (Bruni et al., 2022). Besides, organic fertilizers applied to increase C stocks can change both pH and soil Cu content, leading to supplementary uncertainties (Laurent et al., 2020). Furthermore, together with rainfall and soil moisture changes, climate change is expected to also induce higher temperatures and shorter winters so that a shift in cultures toward the north is expected (Hannah et al., 2013). Therefore, areas with currently low total soil Cu levels may potentially experience a rise in Cu inputs from fungicides, which may subsequently be transported through freshwater systems. Thus, the estimations of LP and AP as computed here can be used to identify regions about to leach or accumulate high amounts of Cu and to anticipate total content modifications that could occur with an eventual change in anthropogenic activities. Indeed, land management changes due to land use changes or regulation changes may affect the use of Cu in agriculture in the future, with potential consequences for Cu leaching.

As a first step, the study conducted here could be used to highlight areas needing regulations to lower Cu input thresholds. Indeed, the changes in the LP (and AP) areas we noticed are not only the reflection of the general runoff change or of the current Cu risk but also underline areas of interest when combining risks linked to soil contamination and climate change. For instance, in eastern Europe, low K_f and high runoff result in Cu LP areas with soils that tend to transfer Cu from soils to 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 LSM and GCM combination and whatever the RCP for at least one period for each RCP examined. In these vineyard regions (Abruzo and Marche regions), annual Cu inputs are high, resulting in Cu accumulation in soil surface horizons. 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, but are often very punctual and local. Hence, to go further, localization of areas with exogenous risks of Cu dissemination must be identified to reinforce the predictions, e.g., by coupling studies of leaching potential, such as the one we conducted here, with erosion risk studies (Panagos et al., 2021) and with outlet characteristics.

4.2 Temporal change in data and scope of the modeling analysis

To reduce intra- and inter-annual variability, the modeling conducted here focused on 5-year means and thus aimed to smooth the seasonal variability of runoff. The K_f we calculated was not a dynamic value since we did not make hypotheses about the temporal change in soil organic carbon or pH. Furthermore, K_f is defined on the assumption that there is equilibrium between the solid and solution phases. This means that the amount of Cu in the solution estimated by this method may be less than that present immediately after Cu application and before equilibrium is reached (McBride et al., 1997). Nevertheless, our results showed good agreement between the four LSM and GCM combinations in terms of their projection of the number of grid cells where both LP and AP are detected, validating the use of their median to perform projections in the absence of in situ validation.

It must be noted that the scope of our predictions has limits that relate to the difficulties in predicting whether rain- and snowfall and runoff will evolve in terms of intensity and frequency. It has already been identified that, during high-load events, much more Cu was transported in the solution than during light events (Imfeld et al., 2020), but alternations between drying and rewetting events may also affect Cu partitioning between phases (Christensen and Christensen, 2003; Han et al., 2001). Also, to gain field reality at the local scale (here, up to 50 km), such as in landscapes or catchments, modeling will require accounting for the time periods of the year with higher rain- and snowfall amounts coinciding with periods of Cu use, for instance, in agriculture and vineyards (Ribolzi et al., 2002; Banas et al., 2010). Indeed, if intense rainfall occurs close to Cu fungicide applications, a larger Cu amount than what is locally computed, taking into account total Cu and K_f, may be exported through runoff (Ma et al., 2006b, a). Thus, local soil Cu budgets require the use of temporal models which account for the regular inputs and outputs of Cu from vegetation and runoff that cannot be accounting for with multi-year means. Finally, the identification of the areas with high risks of soil Cu leaching or accumulation can be viewed as a first step toward the risk change assessment of Cu contamination that would be useful for land management or Cu fertilizer application regulations.