– Modelling the dynamics of the land-sea nutrient transfer over the Mediterranean region–version 1: Model description and evaluation

. Land forcing (water discharge, and nutrient loads) is reported as one of the major sources of uncertainty limiting the capacity of marine biogeochemical models. Runoff from rivers and coastal plains delivers signiﬁcant amounts of nutrients to the Mediterranean Sea from agricultural activities and urban waste water. Several recent studies show that variations in river inputs may play a signiﬁcant role in marine biogeochemical cycles and the planktonic food web throughout the entire basin. The aim 5 of this study is to estimate the water discharge as well as nitrate (NO 3 ) and phosphate (PO 4 ) release into the Mediterranean Sea from basin-wide agriculture and inhabited areas through the implementation of the biogeochemical land-sea nutrient transfer processes within the agro-ecosystem model Lund Potsdam Jena managed Land for the Mediterranean (LPJmL-Med). The representation of the nutrient transfer from land to sea has been introduced into LPJmL-Med by considering the following processes: remineralization, adsorption, nitriﬁcation, denitriﬁcation and phytoplankton dynamics. A compilation of a new input 10 data set of fertilizer, manure and wastewater nutrient content [1961-2005] has been added to the LPJmL-Med forcing data set. The ﬁrst basin-wide LPJmL simulation at 1/12° indicates that the model succeeds in simulating the interannual variability of water discharge for the main rivers in the Mediterranean Sea, especially the Po, Rhone and Ebro Rivers. A very high correlation (R-square values higher than 0.94) is found for these three rivers. Results also show a good consistency between the simulated nutrients concentration (NO 3 and PO 4 ) and available in-situ data. River outﬂows of NO 3 and PO 4 exhibit opposite trends in the 15 Mediterranean Sea. NO 3 showed a more or less continuous increase from the beginning of the 1960s until the present in all three rivers. PO 4 trends are more heterogeneous. There is a strong increase in PO 4 between 1960 and 1980, followed by a decrease in mean annual ﬂuxes from the second half of the 1980s as a consequence of the banning of phosphates in detergents, and the construction of waste water treatment plants in the different countries. Results show that wastewater strongly contributes to the river phosphate ﬂuxes, while both agriculture and wastewater control the nitrogen (mainly as NO 3 ) ﬂuxes from rivers to the 20 Mediterranean Sea.

Biogeochemical cycles are currently undergoing major changes due to human activities, such as intensive agriculture and fossil fuel burning (Gruber and Galloway, 2008;Bouwman et al., 2009;Fowler et al., 2013). The effects of anthropogenic activities on biogeochemical cycles can be either direct (e.g. through an increased supply of nutrients to ecosystems), or 25 indirect (e.g. through the rise in air and water temperature, which in turn alters the rate of key processes in nutrient cycling).
This alteration of the biogeochemical cycles of key nutrients such as carbon, nitrogen and phosphorus may in turn heavily affect the whole structure and functioning of terrestrial and marine ecosystems. For example, the increasing use of nitrogen and phosphorus fertilizers for agricultural production can have devastating effects on both aquatic and marine ecosystems (e.g. Diaz and Rosenberg, 2008;Scavia et al., 2014). The development of cities with increasing populations also affects the release of 30 nutrients through wastewater (Morée et al., 2013), in turn affecting the functioning of aquatic and marine ecosystems. However, nutrient retention along the aquatic continuum can be substantial and have a buffer effect regarding the human-induced increase in nutrient supply to terrestrial and aquatic systems, although it does not balance nutrient release at global scale (Beusen et al., 2016;Lassaletta et al., 2012;Billen et al., 2007). Retention processes also modify the ratio at which nutrients are transported downstream, as well as their chemical form, which can have important implications for biogeochemical cycles from local to 35 global scale Peñuelas et al., 2013). Assessing the amount of nutrients retained in soils and rivers is thus crucial to understand the impact of human-induced release of nutrients on aquatic and marine ecosystems.
Changes in agriculture in industrialized countries over the last 50 years have strongly affected the severity of the impact of anthropogenic activities on global biogeochemical cycles, mainly with regard to the nitrogen and phosphorus cycles Fowler et al., 2013). Switching from traditional to modern agriculture has resulted in a geographical separation of 40 livestock and crop production, inducing an increased use of inorganic fertilizers to substitute for manure (Thieu et al., 2011).
In areas with intensive livestock production, the content of nitrogen and phosphorus in soils has strongly increased, inducing major extensive leaching of these nutrients in soils and the contamination of freshwater (Thieu et al., 2011;Bouwman et al., 2013;Leip et al., 2015). In cropping systems, a large part of nutrients from manure and inorganic fertilizers is not assimilated by plants and is thus transported through leaching and runoff along the aquatic continuum (Coskun et al., 2017). Eutrophication 45 in coastal waters due to the increased inputs of nitrogen and phosphorus from groundwaters and rivers may indeed induce hypoxia or anoxia, and can lead to changes in the whole trophic web (Diaz and Rosenberg, 2008;Rabalais et al., 2010;Moller et al., 2014).
The Mediterranean region includes countries with a wide range of socio-economic and development backgrounds. A specificity of this region is the diversity of the agricultural practices in the different countries, characterized by differences in crop 50 types, fertilization and irrigation methods, as well as livestock and manure quality (Rost et al., 2008;Potter et al., 2010;Fader et al., 2015). The Mediterranean region also includes a wide variety of rivers, which differ by their flow rate, their seasonality, their degree of exploitation for irrigation, and wastewater treatment policy. The region is subject to global change pressures such as increasing population and gross domestic product (GDP), urbanization, and climate change (Schröter et al., 2005;Diffenbaugh and Giorgi, 2012;Zdruli, 2014). Its agricultural production will probably face higher vulnerability both with re-gard to ensuring food security (several countries already requiring extensive food imports), and for exporting products that are essential to the national economies of many countries (Hervieu, 2006;Verner, 2012). Thus, it is important to understand how the effects of urbanization, socio-economic and agricultural development in Mediterranean landscapes translate into land-sea nutrient transfer dynamics.
The most important macro-nutrients for the marine ecosystem are nitrate and phosphate. Nitrogen (N) and phosphorous 60 (P) occur in rivers in various forms, i.e. dissolved, particulate, organic or inorganic. Inorganic form is becoming increasingly important due to the growing anthropogenic contribution. Although N is generally considered to limit primary productivity in most of the world's oceans, previous studies have suggested that the Mediterranean Sea may be an exception, with a strong P limitation (Krom et al., 1991). The degree of P limitation increases from west to east across the entire basin (e.g. Krom et al., 2003). For these reasons, and because of the greater availability of data on NO 3 and PO 4 as compared to the other nutrient important characteristic of the LPJmL model is to be dynamic, as the properties of crops can change temporally and are locally dependent on external forcing. For example, agricultural practices such as sowing date depend on climatic conditions (which impacts the cultivar choice), and crop management as irrigation practices and residue management can be modified (Bondeau 100 et al., 2007).
LPJmL basic version uses 13 CFTs (11 annual arable crops and two managed grass types), with specific parameterizations (Bondeau et al., 2007). Additional CFTs can be added for specific applications (Lapola et al., 2009;Jans et al., 2020). For Mediterranean studies, Fader et al. (2015) recently included 10 additional crop types in LPJmL, mainly perennial crops, which are important in the Mediterranean and are often major consumers of irrigation water. With this LPJmL-Med version, the model 105 can estimate the impacts of different demographic projections and irrigation management systems on production and irrigation water consumption under climate change for the entire Mediterranean basin (Fader et al., 2016).
So far, LPJmL has not explicitly accounted for the release of nutrients (NO 3 and PO 4 ) to the sea from agricultural activities and urban waste water (via runoff from rivers). In the present study the representation of the nutrient transfer from land to sea is introduced into LPJmL-Med by considering the following processes: adsorption, remineralization, nitrification, denitrification, 110 and phytoplankton dynamics (Fig. 1). A compilation of a new input data set of fertilizer, manure and wastewater nutrients content, at an annual time step for the period  has been established at high spatial resolution grid (1/12°) and over watersheds of rivers flowing into the Mediterranean Sea (Fig. 2). All implemented processes and the corresponding references are described in detail in the following sections. 115 We describe here all new model features in full detail, i.e the representation of the nutrient transfer from land to sea, by considering the following processes: adsorption, remineralization, nitrification, denitrification, phytoplankton dynamics, and by the implementation of new external input of nutrients. Therefore, we start with general equations describing the nutrient dynamic in soil and in water (cf. below). All parameters and variables used in the model are summarized in Tables A1, A2 and A3 in the Appendix.

General equations
State variables used to describe C, N and P cycles in the LPJmL-Med model are the following: -Residues content in the soil layer l (RES C,l , RES N,l and RES P,l in g).
-Decomposers content in the soil layer l (DEC C,l , in g ), only represented through a carbon content, their N and P contents being deduced from fixed C:N and C:P ratios.
-Nitrate, ammonium and phosphate content in layer l of the soil (NO3 l , NH4 l , PO4 l in g).
-Detrital organic matter concentration in water (RES C,water in gC·m −3 ).

Soil residues 130
The soil residue pool consists of the organic matter coming either from the dead biomass entering the litter pool or from the applied manure. Its mineralization will produce inorganic nutrients. The residue pool is updated for each soil layer, at different time steps depending on the source of the organic matter. Manure is applied in agricultural fields twice a year. For natural PFTs, the litter pool is incremented monthly with the fine root turnover, daily (for evergreen PFTs) or at leaf fall (for deciduous PFTs) with the leaf turnover, yearly with the wood biomass of killed trees. For CFTs, the litter pool is incremented at harvest 135 with fractions of the different organs depending on the farming practices. Residues in the soil layers below the surface layer are fed by roots only. Part of the residues is consumed by soil decomposers.
The dynamics of carbon, nitrogen and phosphorus (X either stands for N or P) contents of soil residues is described by the following equations for the above soil layer (l = 0): and for the deeper soil layers, in contrast to residues added to the first soil layer, those added in underneath soil layers (1 ≤ l ≤ 4) only come from roots : where VEG C,i is the C biomass of the plant part i (i = root, sapwood, leaf, storage organ) in the top soil layer and ROOT C,l , 150 the C biomass of roots in layer l (1 ≤ l ≤ 4). See section 2.4.1 for more details. FRAC res,i is the daily fraction of plant part i entering the residue compartment of the above soil layer due to senescence, mortality, or harvest (section 2.4.1). FRAC res,root,l (l > 0) is the daily fraction of root entering the residue compartment of soil layer l due to senescence, mortality, or harvest.
RATIO C:X,i is the C:X ratio for the plant part i. MAN X is the content in element X of the daily applied manure (section 2.3.2). G dec,l is the daily amount of C residues consumed by soil decomposers in layer l (section 2.4.1).

Soil decomposers dynamics
Changes in the C biomass of soil decomposers in each soil layer l (0 ≤ l ≤ 4) is assessed as follows: where m dec is the daily mortality rate of decomposers (cf. Table. A1).

Water residues dynamics
Detrital organic matter consumed by decomposers are substracted from the residue pool in the water: where PHYTO C,water is the Phytoplankton C biomass (cf. equation 11), m phyto is the daily mortality of phytoplankton 165 and water decomposers, including grazing, sed is the fraction of residues trapped daily into sediments, G dec,water is the daily amount of C residues consumed by water decomposers (section 2.4.2).

Water decomposers dynamics
Changes in the C biomass of water decomposers (DEC C,water ) is assessed as follows: ∂DEC C,water ∂t = G dec,water − m phyto · DEC C,water Nutrient dynamics in soil (daily time step): The nutrients (i.e. NO 3 , NH 4 and PO 4 ) dynamics in soil can be represented by the following equations: LEACH NO3,l , LEACH NH4,l , and LEACH PO4,l : leaching, i.e. transfer of nutrient in the different soil layers = inflow -outflow of 185 nutrient through water percolation and surface/subsurface runoff (Sect 3.2).

Phytoplankton dynamics
Phytoplankton C biomass (PHYTO C in gC) is calculated in equation 11 as: G phyto is the daily growth rate of Phytoplankton (cf. section 2.4.2) and m phyto is the daily mortality of phytoplankton and 190 water decomposers, including grazing.
Nutrient dynamics in water (daily time step): The nutrients (i.e. NO 3 , NH 4 and PO 4 ) dynamics in water can be represented by the following equations: ∂NO3 water ∂t = SEW NO3,water − UPT NO3,water + NITR NH4,water − DEN NO3,water + TRANS NO3,water Note: (i) Nitrogen-fixation is implicitly represented by considering the N uptake from N-fixers (the legumes soybeans and pulses) comes from N-fixers only, resulting in UPT NO3 = UPT NH4 = 0 for those CFTs. Therefore, there is no specific Nitrogenfixation element within the global equations (ii) Although NO x deposition has clearly increased in the last decades, its values in the Mediterranean remain low (apart in the Po region) when compared to the load from agriculture. It is not included in the 210 equations.

Implementation of external inputs of nutrients
In this section, we will focus on processes leading to external inputs of nutrients to soils and rivers, either associated with agricultural practices (i.e., application of manure and inorganic fertilizers to soils) or with household wastewater release.

215
N and P fertilizers are applied to the first soil layer and depend on the country, the crop type and the cultivated area. N and P fertilizers are applied twice or thrice each year depending on the crop fertilization calendar ( Table 1). The IFADATA database (IFA, 2016) provides N and P fertilizer data, we consider that half of the inorganic N inputs from fertilizers enters the nitrate pool in the first soil layer, and the remaining half feeds the ammonium pool. On fertilizer days, we estimate the amount of , NH4, and PO4fertilizers applied to each stand on the first soil layer (l=0), though equation 15 as in Potter et al. (2010): FERT N,IFA and FERT P,IFA are the annual inputs of N and P fertilizers obtained from the IFADATA database IFA (2016) for each country. AREA is the area of the stand, while AREATOT is the total cultivated area in the country, both obtained from the land use input data of the LPJmL-Med model (cf. Fig. 3). FRAC f ert corresponds to the fraction of fertilizers applied each fertilizer day (1/2, 1/3 or 0 depending of crop type, 1).

Manure application
As for inorganic fertilizers, manure is applied to the first soil layer and depends on the country, the crop type and the cultivated area. As it is composed of organic matter, manure is added to the litter (see Fig. 1). To estimate the amount of N manure LPJmL-Med model (Fig. 4). In contrast to inorganic fertilizers, we consider that manure is applied twice each year for all crop types except pastures where it is supposed to be spread) uniformly over the grazing period 6 months long in the northern hemisphere, all year long in the southern hemisphere, as described in the crop fertilization calendar (Table 1). As data for P manure application are not available from the FAOSTAT database, we estimate it using a constant country-specific

Wastewater release
Wastewater is released into rivers and reservoirs ( Fig. 1, Fig. 5). We distinguish the emissions by urban and rural populations, 245 as they can have an uneven access to wastewater treatment plants (Van Drecht et al., 2003). The daily point source of N and P from sewer systems (SEW X , X =N, P) is estimated following equations 18 from Van Drecht et al. (2003: WW X,urban and WW X,rural (X =N, P) correspond to the annual per capita X release from wastewater in urban and rural areas. Urban and rural population size in each grid cell (POP urban and POP rural , respectively) are provided each ten years The N content of wastewater is determined only by human emissions (E N,hum ), while P release is determined by both human emissions (E P,hum ) and emissions from laundry and dishwasher detergent (E P,Ldet and E P,Ddet , respectively) as shown in 255 equation 20 (Van Drecht et al., 2003: FRAC sew,i WW P,urban = E P,hum + E P,Ldet + E P,Ddet FRAC sew,urban · FRAC sew,urban · (1 − REM P,1 ) WW P,rural = E P,hum + E P,Ldet + E P,Ddet FRAC sew,rural · FRAC sew,rural · (1 − REM P,1 ) 260 FRAC sew,urban and FRAC sew,rural in equations 20 correspond to the fraction of the urban and rural population that is connected to public sewerage system. Four categories of treatment are considered: no treatment (i = 0), primary/mechanistic treatment (i = 1), secondary/biological treatment (i = 2), and tertiary/advanced treatment (i = 3). The fraction of the total population connected to public sewerage system with treatment i is determined by FRAC sew,i . REM N,i and REM P,i correspond to the fractions of N and P removed from wastewater with treatment i (cf . Table A1). When no data is available concerning the 265 access of population to the different types of sewerage systems, we consider that population can only access primary treatment.
The per capita human emissions of N (E N,hum ) and P (E P,hum ) are estimated from the national per capita gross domestic product based on purchasing power parity (GDP P P P ) as shown in equations 21 (Van Drecht et al., 2003: ww where RATIO N :P,ww is the mean N:P ratio of municipal wastewater. The per capita emission of P from laundry and dishwasher detergents (E P,Ldet and E P,Ddet , respectively in equations 22) are estimated from socioeconomic indices (Van Drecht et al., 2009): CONT P,Ldet and CONT P,Ddet are the P content of laundry and dishwasher detergents. FRAC P f ree,Ldet is the fraction of P-free laundry detergents. The fraction of P-free dishwasher detergents is assumed to be negligible (Van Drecht et al., 2009). GDP mer corresponds to the national per capita gross domestic product based on market exchange rate. CONS hh,Ddet corresponds to the mean consumption of dishwasher detergents for each household that owns a dishwasher. PPHH is the 280 average number of persons by household.
The fraction of P-free laundry detergent is estimated in equation 23 as in Van Drecht et al. (2009), and depends on the period concerned: -In year yr between 1960 and 1999: -In 2000 and afterwards: FRAC P f ree,Ldet (yr) = min GDP mer 33000 , 1

Transformation and retention processes 290
In this section, we will focus on processes leading either to a change in the form of nutrients (i.e. primary production, residue decomposition and nitrification) or to a loss of nutrients from soils and rivers (i.e. sedimentation, adsorption and denitrification).

Primary production and residue decomposition in terrestrial systems
Both natural and cultivated plants consume nitrate, ammonium and phosphate in soils. Unlike more complex modelling of the nutrient uptake (e.g. Von  for N), we consider here a simple formulation, similar for N and P, in order to estimate 295 their daily uptake by plants. As the growth of all plant parts (i.e. leaves, roots, sapwood and storage parts) is already simulated in the LPJmL-Med model in terms of C, we use the specific stoichiometry of the different plant parts to infer the corresponding daily uptake of nitrate, ammonium and phosphate (UPT N O3 , UPT N H4 and UPT P O4 , respectively, cf equations 25). In the LPJmL-Med model, the soil is divided in 5 layers (numbered from 0 to 4) in order to represent the vertical heterogeneity in water and carbon content of soils. Nutrient uptake by plants in the soil layer l is then assumed to depend on the fraction of total 300 roots of this soil layer (FRAC root,l ):

305
UPT C,i is the daily C uptake for the growth of plant part i (i = root, sapwood, leaf, storage) derived from the LPJmL-Med model. RATIO C:N,i and RATIO C:P,i correspond to the C:N and C:P ratios of plant part i (Mooshammer et al., 2014). Note that for crops associated with N-fixers (i.e. soybean and pulses), all the N uptake is assumed to come from symbiotic N-fixers, Part of the plant biomass is added to the soil residue pool at senescence. This can also be the case at harvest depending on 310 the crop residue management. Residues added to the first soil layer (l = 0) correspond to the roots in that layer, as well as to dead organic material from the above-ground parts of the plants (equations 26): VEG X,i is the X (X = C, N, P ) biomass of the plant part i (i = root, sapwood, leaf, storage organ). FRAC res,i is the daily fraction of plant part i entering the residue compartment of the above soil layer due to senescence, mortality, or harvest. SAP C , 320 LEAF C and STOR C are the total amount of C in roots, sapwood, leaves and storage parts of the plant, respectively. ROOT C is the turnover fraction of C in root, FRAC res,i corresponds to the fraction of plant part i entering the residue compartment due to senescence or harvest. In the case of non-permanent crops, FRAC res,i depends on residue management practices.
Part of the residues is consumed by soil decomposers. We assume that residue consumption in soils follows a donorcontrolled functional response (Zheng et al., 1997), and thus is independent of decomposer biomass. As soil decomposers 325 have a different stoichiometry from residues, their growth will be limited by the scarcest nutrient (either N, P or C). In order to maintain homeostasy (e.g. Daufresne and Loreau, 2001), soil decomposers can immobilize part of the inorganic nutrients available in the soil (equation 27): G dec,l = min VMAX res,l S(24, 14) RES C,l , G dec,l is the daily amount of C residues consumed by soil decomposers in layer l (1 ≤ l ≤ 4). VMAX res,l corresponds to the maximum fraction of C residues consumed daily by soil decomposers (Schjonning et al., 2004). S(T opt , γ) defines a bell shaped temperature dependence of soil decomposition (see Billen et al. 1994), as described in equation 28. IMMO X is the maximum fraction of inorganic nutrient X (X = N, P) immobilized daily by soil decomposers to meet their nutrient requirements. RATIO C:N,dec and RATIO C:P,dec are the C:N and C:P ratios of soil decomposers, respectively.
The bell shaped temperature-dependence of biological processes is determined as follows: where T opt is the optimal temperature for the process considered, γ is the sigmoid range of T, and T env is the environment temperature (i.e. soil or water temperature depending on the process considered).

340
Finally, the remineralization REMIN X,l , (X = NH 4 and PO 4 ) following the decomposition of organic N and P is determined as follows: where m dec is the daily mortality of decomposers.

Primary production, sedimentation and residue decomposition in rivers and reservoirs
Phytoplankton consumes nitrate, ammonium and phosphates in rivers and reservoirs. We consider that phytoplankton consumes nitrate and ammonium in the same proportion as they are available in water. We assume that Liebig's law of the minimum (Von Liebig, 1942) governs the growth of phytoplankton, i.e. that its daily growth rate (G phyto ) is limited by either N or P as shown in equation 30: µ phyto is the maximum daily growth of phytoplankton (in d −1 ). K N,phyto in gN/m 3 and K P,phyto in gP/m 3 correspond to the half-saturation constants for phytoplankton growth regarding N and P, respectively. S(T opt, γ) is already described in equation 28. V water is the volume of water in the river or reservoir considered (in m 3 ). NO3 water , NH4 water and PO4 water are the concentrations of nitrate, ammonium and phosphate, respectively, in water.

355
Part of the remaining residues are consumed by water decomposers. We assume that organic matter consumption in water follows a donor-controlled functional response, and thus is independent of decomposer biomass (see equation 31): G dec,water = VMAX res,water · S(24, 14) · RES water (31) G dec,water in equation 31 is the daily consumption of carbon residues by water decomposers. VMAX res,water corresponds to the maximum fraction of residues consumed daily by water decomposers.

360
The uptake in inorganic nutrient content of water by phytoplankton is calculated by the following equations 32: 365 where RATIO C:N,phyto and RATIO C:P,phyto are the C:N and C:P mean ratios of phytoplankton cells. Note that we consider that water decomposers and phytoplankton have the same characteristics regarding their mortality rate and stoichiometry.
The remineralization in water REMIN X,water , (X = NH 4 and PO 4 ) following the phytoplankton water decomposers is determined as follows: Where: a and b are fixed parameters (cf. table A1). m l is the mass of soil in layer l. VOL water,l corresponds to the volume of water in the soil.

380
In rivers, we consider that the amount of NH 4 adsorbed on particles is negligible. The daily amount of PO 4 adsorbed in rivers and reservoirs (ADS P O4,water ) follows a Michaelis-Menten dynamics and depends on the concentration of suspended particles in the water (C susp,water ), as described in Billen et al. (2007) (see equation 35): ADS P O4,water = C susp,water · VOL water · VMAX P O4,ads · [PO4] water [PO4] water + K water,P O4 · VOL water VMAX P O4,ads corresponds to the maximal adsorption rate of P O 4 on suspended particles in rivers and reservoirs. Kads water,P O4 385 is the half saturation constant for P O 4 adsorption in rivers and reservoirs.
As land use impacts the stability of soils and the amount of organic material released in aquatic systems, it also affects C susp,water the concentration of suspended particles in rivers and reservoirs (Billen et al., 2007) as shown in equation 36: C susp,i is the concentration of suspended solids in the water released from system i (i = forests, grasslands, crops, urban areas),

390
SURF i corresponds to the surface of system i in the cell, whose total surface is SURF cell . The amounts of P O 4 lost through adsorption in soils and water is then substracted from their respective P O 4 pool.

Nitrification
Nitrification is an important process in both terrestrial and aquatic nitrogen cycles, corresponding to the transformation of NH 4 into NO 3 . In each soil layer l(l = 0 − 4), the daily amount of NH 4 transformed through nitrification (NITR l in gN) depends on 395 soil moisture, temperature and pH Parton et al., 1996): NITR l = NH4 l · VMAX soil,nitr · RESP nitr,moist · S(18.79, 7.44) · RESP nitr,pH VMAX soil,nitr is the maximum fraction of NH 4 nitrified daily. RESP nitr,moist and RESP nitr,pH correspond to the response functions of soil nitrification regarding soil moisture and pH, respectively. S(18.79, 7.44) defines a bell shaped temperature dependence of soil nitrification, as described in equation 28 Parton et al., 1996).
FRAC soil,water corresponds to the ratio of soil moisture to soil maximum moisture. a, b, c and d (0.60, 1.27, 0.0012, and 2.84 respectively) are fixed parameters for sandy and medium soils (Von Parton et al., 1996).

405
The response function of nitrification regarding soil pH (pH soil ) is based on Parton et al. (1996) (see equation 39): In rivers and reservoirs, nitrification is modeled through a Michealis-Menten equation, as in the Riverstrahler model (Billen et al., 1994), equation 40: NITR water = VMAX water,nitr · S(24, 6) · NH4 water NH4 water + K water,nitr NITR water corresponds to the daily amount of NH 4 nitrified in water. VMAX water,nitr is the maximum nitrification rate in water. K water,nitr is the half-saturation constant for water nitrification. S(24, 6) defines a bell shaped temperature dependence of nitrification (Billen et al., 1994), as described in equation 28.
The amount of N nitrified in soils and water is then added to their respective NO 3 pools and substracted from the NH 4 pools.

415
Denitrification is a key process occuring along the land to sea N transfer, as it allows the release of N into the atmosphere in the form of dinitrogen. It occurs in hypoxic or anoxic areas, both in soils and aquatic systems. As oxygen availability is not explicitly represented in our model, we use proxies to determine the strength of denitrification in soils and water.
In each soil layer l(l = 0 − 4), the daily amount of N lost through denitrification (DEN l ) depends on soil moisture, temperature and organic content (Von Parton et al., 1996), equation 41: where: VMAX den,l is the maximum daily denitrification in soil. RESP den,moist and RESP den,temp are the response functions of soil denitrification regarding soil moisture and soil temperature, respectively. RES C,l is the mass of carbon in residues of the soil layer l.

425
In soils, we consider that the oxygen content of a soil layer is inversely related to soil moisture (Von Parton et al., 1996): RESP den,moist = 6.664096 · 10 −10 · exp (21.12912 · FRAC soil,water ) The response function of soil denitrification regarding soil temperature (T soil ) is assumed to be ascending up to a temperature threshold of 45.9 • C, and zero beyond (Von Parton et al., 1996), as suggested by equation 43: In rivers and reservoirs, we consider that the oxygen content of water is positively correlated to the hydraulic load (LOAD water in m/d), thus the strength of denitrification varies inversely with the hydraulic load (Wollheim et al., 2008): DEN NO3,water = NO3 water · VMAX den,water lost through denitrification in rivers and reservoirs. S(24, 14) defines a bell shaped temperature dependence of denitrification (Billen et al., 1994), as described in equation 28. VIT water,den is the uptake rate for denitrification in rivers and reservoirs. The amount of N denitrified in soils and water is then substracted from their respective NO 3 pools.

Rivers, dams, and lakes
LPJmL-Med divides the soil column into five hydrological active layers of 0.2, 0.3, 0.5, 1, and 1 m thickness (Sibyll Schaphoff et al., 2013). Water holding capacity and hydraulic conductivity are derived for each grid cell using soil texture from the Harmonized World Soil Database (Nachtergaele et al., 2014) and relationships between texture and hydraulic properties from Cosby et al. (1984). Water content in soil layers is altered by infiltrating rainfall, gravity (percolation), and the plant water 445 uptake. The infiltration rate of rain and irrigation water into the soil depends on the current soil water content of the first layer. The surplus water that does not infiltrate is assumed to generate surface run-off. The lateral exchange of water discharge between grid cells through the river is computed via the river routine module implemented in LPJmL-Med model. The transport of water in the river channel is approximated by a cascade of linear reservoirs Rost et al., 2008).
The human influence on the hydrological cycle is explicitly represented in LPJmL-Med model by accounting for irrigation, 450 water consumption, water abstraction, as well as an implementation of dams and reservoirs. Reservoirs are filled daily with discharge from upstream locations and with local precipitation. Cells receive water from the reservoirs when the following conditions are met: (i) cells have a lower altitude than the cell containing the reservoir, and (ii) they are situated along the main river downstream or at maximum five cells upstream. Hence, a cell can receive water from multiple reservoirs. Dams built for irrigation are assumed to release their water proportionally to gross irrigation water demand downstream. Dams built 455 for other purposes (e.g flood control, hydro-power) are assumed to release a constant water volume throughout the year. The actual release from a reservoir is simulated to depend on its storage capacity relative to its inflow. If an irrigation purpose is defined for the reservoir, part of the outflow is diverted to irrigated lands downstream. Both surface and subsurface run-off are simulated to accumulate to river discharge (for more detail see Schaphoff et al., 2018).

Nutrient leaching and transport 460
Nutrients are assumed to be fully dissolved in water and move with lateral runoff, surface runoff, and percolation water. The first step to calculate the quantity of nutrients transported with the water from a soil layer, is to compute the concentration of nutrients in the mobile water. This concentration is then multiplied by the volume of runoff or percolation water between soil layers or into the aquifer. The quantity of nutrients leached depends on climate conditions (e.g precipitation), soil conditions and the intensity of soil management (e.g fertilization, plant cover, soil treatment). Nutrient movement with water fluxes is simulated as in SWAT (Neitsch et al., 2002(Neitsch et al., , 2005, and in Von  in LPJmL5.0 version. The concentration of nutrients (Nut) in the mobile water (in kg m −3 ) is calculated from equation 45: w mobile is the amount of mobile water in the layer (mm), θ = 0.4 is the fraction of porosity from which anions are excluded (the same as in Von , 0.5 in Neitsch et al. 2002, and SAT l is the saturated water content of the soil layer (mm).

470
The mobile water w mobile,l in the layer l is the quantity of water lost by surface runoff, lateral flow, and percolation as shown in equation 46: where Q surf is the surface runoff (only in the top soil layer; mm), Q lat,l is the water discharged from the layer by lateral flow (mm), and w perc,l is the amount of water percolating to the underlying soil layer on a given day (mm).

475
Finally, the flux of nitrate that is removed through surface runoff Fnut surf and lateral flow FNut lat,l is calculated in equation 47 as: Fnut lat,l=0 = β N ut · NUT mobile,l=0 · Q lat,l=0 Fnut lat,l = β N ut · NUT mobile,l · Q lat,l

480
In deep soil layers, β N ut which is the nutrients percolation coefficient, it controls the amount of NUT removed from the surface layer in runoff relative to the amount removed via percolation Neitsch et al. (2002). The value for β N ut can range from 0.01 to 1.0. We choose for β N ut a value of 0.4 (similar to Von Bloh et al. (2018)).
Nutrients flux moved to the lower soil layer with percolation Fnut perc,l is calculated in equation 48 as: We start by spin-up simulation, where water pool and carbon are initialized to zero during 5000 years to bring natural vegetation patterns, and carbon stocks into dynamic equilibrium. We cyclically repeat 30 years of climate data input with constant concentrations of atmospheric carbon dioxide (at 278 ppm). Then a second phase of spin-up (during 500 years) during which land use is introduced in the year 1700 from which it is updated annually according to the historic land use data-set (see Fader et al., 2015, 2010, andSect. 4.1). And finally we run a hind cast simulation  with the new external input of nutrients 495 described in the following section.

Model input
The model domain encompasses the catchments of rivers ending in the Mediterranean Sea (see Fig. 2). We must note that, since the assessment of some of the simulated variables requires data available at the administrative level, we run the model over the entire area of the Mediterranean countries, i.e. also for catchments not ending in the Mediterranean Sea. The required input 500 data for LPJmL-Med are: (i) gridded climate variables (temperature, precipitation, and radiation); (ii) atmospheric CO 2 concentrations; (iii) gridded soil texture as described in Schaphoff et al. (2018); (iv) gridded land use and crop distribution dataset; and (v) for the present application we have implemented three new inputs quantifying the nutrients load from agricultural management and wastewater release (see below).

505
The model runs on a daily time step. For the present study, we used monthly climate inputs data (i.e. precipitation and temperature data) from the regional climate model CNRM-ALADIN (Herrmann et al., 2011;Farda et al., 2010;Déqué and Somot, 2008) in the framework of Med-CORDEX programme (https://www.medcordex.eu/). Here we use the high-resolution version at 12 km (converted to a resolution of 5 by 5 arc-minutes) for the area covered by the domain of the CNRM-ALADIN model (i.e. between 27°N and 57°N in latitudes and between 10°W and 40°E in longitudes). Monthly temperature is linearly interpo-510 lated, while a generator based on monthly precipitations and monthly number of wet days provide daily precipitation values (Gerten et al., 2004). To complete the missing ares (i.e. southern 27°N) we have used the CNRM-CM5 global model (Voldoire et al., 2013); from the ESGF data set, https://esgf-data.dkrz.de/projects/esgf-dkrz/). Incident solar radiation is internally computed from the cloudiness data of the Climate Research Unit's (CRU) time series (TS) 3.1 data (Mitchell and Jones, 2005).
LPJmL works with soil hydrological parameters that depend on soil texture. We use the 13 USDA soil texture classes that are 515 provided by the Harmonized World Soil Database v 1.2 (FAO, 2012) at 30 arc-second spatial resolution. The original data have been agregated at 5 arc-minute resolution.

Land-use and river network data
The land use data for the crops in LPJmL-Med had been compiled from different sources (Portmann et al., 2010;Monfreda et al., 2008;Klein Goldewijk et al., 2011), as explained in Fader et al. (2015Fader et al. ( , 2010. Decadal cropland data from HYDE were interpolated to derive annual values and then used for extrapolating the land use patterns of ∼2000 to the past, until 1700.
Historical irrigation fractions were determined as explained in Fader et al. (2010). Further information is given by Fader et al. (2010Fader et al. ( , 2015. Runoff is simulated by the model for each grid cell, and a drainage direction map gives the transport directions of the flows towards the different rivers. The river-routine scheme is derived from the Hydrosheds database (Lehner et al., 2008) 525 and is updated at 5 arc minute resolution (Siderius et al., 2018). The GRanD database (Lehner et al., 2011) provides detailed information on water reservoirs that includes information on storage capacity, total area, and main purpose. Furthermore, information on natural lakes is obtained from Lehner and Döll (2004).

Manure inputs
The Food and Agriculture Organization provides country-specific yearly data on N manure applied to cultivated soils or left 540 on pasture (FAO, 2016a) since 1961. As for fertilizers, adjustments are done when country boundaries have changed (cf.
Appendix B). N:P ratios in manure derived from Potter et al. (2010Potter et al. ( , 2011a are used for estimating P manure applied to cultivated soils or left on pasture. Unlike fertilizers providing mineral N and P nutrients, manure provides organic nutrients that are not immediately available for the plants, therefore manure is applied earlier in the year than the synthetic fertilizers. Our crop fertilization calendar considers two annual application dates for all crop categories (cf. Table 1). For pasture, manure 545 is equally distributed between May 1st and September 30 th . Maps illustrating the temporal dynamic of manure inputs over the pasture and crop areas in the Mediterranean countries is shown in Fig. 4. The equations used for calculating the annual nitrogen and phosphorus content of Manure are described in Section 2.3.2.

Wastewater inputs
Wastewater nutrient content depends on population connected to public sewerage system and on country-specific Gross Domestic Product (GDP). High GDP generally implies a high level of households owning washing machines and dishwashers, therefore a high consumption of detergents (containing more or less pollutants depending on regulations), but it is also associated with a more widespread access to sewerage systems with wastewater treatment. Thus, depending on the demography dynamics and on the GDP change, the wastewater inputs can increase or decrease. The equations and data used for calculating the annual nitrogen and phosphorus content of wastewater are described in Section 2.3.3. Maps illustrating the temporal 555 dynamics of wastewater inputs over the Mediterranean countries is shown in Fig. 5.

Model outputs
As a function of agricultural management and climatic conditions, LPJmL-Med simulates, spatially explicitly and at a daily to yearly temporal resolution, growing periods (sowing and harvest dates), gross and net primary productivity, carbon stock in plants, litter, and soil, agricultural production, as well as a number of hydrological variables, such as soil evaporation, 560 infiltration, percolation, water stress, irrigation water requirements, runoff, and river discharge. In the frame of the present study we have added a new output to the LPJmL-Med model: nutrient concentrations (PO 4 and NH 4 , and NO 3 ) in surface and deep soil, in lake water, and in water discharge, as well as phytoplankton content in lake and in water discharge.

Results and discussion
In the following section, the capacity of LPJmL-Med for simulating the water discharge and nutrient concentrations was 565 compared to published in-situ data, with a focus on the main rivers of the Mediterranean Sea. In general, data coverage of the northern-European rivers is good, whereas data on southern-Mediterranean rivers is poorer. Particularly, long time series of data are missing for southern-Mediterranean rivers.

Water discharges from the main rivers of the Mediterranean Sea
Water discharge is the main factor controlling nutrient transfer from land to sea by rivers. We therefore start our assessment of 570 LPJmL-Med with an evaluation of the rivers water discharge to the Mediterranean. Moreover, the freshwater discharges into the sea can have an important influence on the functioning of marine ecosystems through their control of the ocean dynamics in the Mediterranean Sea (e.g. Skliris et al., 2007). For all these reasons, it is essential to accurately simulate river discharge. This is particularly true in the Mediterranean context, where water resources represent an important economic value and political barriers probably still inhibit greater transparency regarding river data. We chose for the comparison the period between 1920 575 and 1985 because data for the main rivers (i.e., Rhone, Po, Ebro, but not Nile) are available for this period, and to limit the impact of damming and anthropogenic water use. Figure 6 shows the water discharge time series for the Rhone, Po, Ebro and Nile rivers against historical in-situ data between 1920 and 1985. The LPJmL-Med simulates relatively well the interannual variability of the main rivers in the Mediterranean Sea especially the Po, Rhone and Ebro. A very high correlation was found between model outputs and in-situ data (Fig. 7), 580 i.e. the R-square value is higher than 0.95 for the three rivers, showing that the interannual variability is clearly well simulated by the LPJmL-Med model ( Fig. 6 and Fig. 7). Almost all of them presented a significant discharge decrease as a result of the combined effect of ongoing climate change and enhanced anthropogenic water use (Ludwig et al., 2009(Ludwig et al., , 2003. The Ebro river show a net decrease of water discharge between 1960 and 1990 (from 28 to 9 km 3 .y −1 ) due to the higher frequency of dry periods observed in the Ebro River basin during the same period (Valencia et al., 2010). However, the general trend for the Po 585 and the Rhone rivers is to remain about constant for the same period. This distinctive behaviour pattern may be explained by the non-Mediterranean climate in the upper northern parts of their basins (Lutz et al., 2016;Ludwig et al., 2009). The succession of dry and wet periods that can be seen on in-situ data for the three rivers is well reproduced by the model, indicating that the seasonal variability of water discharge is simulated by the LPJmL-Med. However, as shown in Figure 8, the amplitude of the seasonal cycle is larger in the model because we do not represent explicitly the role of dams in the regulation of water 590 flows (i.e., release and retention of water). The construction of dams in the context of water resources management plays an important role matching human needs with the hydrological regime, i.e., to ensure an adequate supply of water by storing water in times of surplus and releasing it in times of scarcity (Fig. 8). This has been well studied for the Ebro catchment (Radinger et al., 2018). It is difficult to represent explicitly this water regulation by dams because very few data are available and those regulations are sporadic and hard to characterize, particularly for long periods.

595
It is also very difficult to evaluate and simulate the discharge from the Nile, which loses most of its water through infiltration in swamps, river evaporation and anthropogenic water use (Nixon, 2003). Since 1965, the construction of the Aswan High Dam has had a major impact on the water discharge. The flow rate of the Nile river is between 6 km 3 yr. −1 (ElElla, 1993) and 15 km 3 .yr −1 (Nixon, 2003). In another study, Skliris et al. (2007) have estimated the flow of the natural long-term discharge at Aswan at about 83 km 3 .yr −1 . However, the Aswan dam is still very far from the Mediterranean Sea and the anthropogenic 600 water use in this area is very low compared to the northern part of Egypt in the Nile delta. In this context, LPJmL-Med simulates relatively well the water discharge of Nile with an average value of 17.88 km 3 .yr. −1 (average value between 1965 and 1985), despite being relatively higher than the previous estimation from literature (ElElla, 1993;Nixon, 2003), but still far lower than the estimation of Skliris et al. (2007) of about 83 km 3 .yr −1 at Aswan.  (Fig.9a), Po (Fig. 9b) and Ebro (Fig.9c) rivers increased steadily from the beginning of the 1960s up to the 1990s. The model simulates a rapid increase in NO 3 in recent years after 1994, with a greater interannual variability for the three rivers. The same increase in NO 3 was presented in Ludwig et al. (2009Ludwig et al. ( , 2003 610 combining in-situ data and the NEWS-DIN model of Dumont et al. (2005), and also from the MOOSE (Mediterranean ocean observing system on environment) observatory system (Raimbault and Lagadec, 2012) available between 1980 and 2004 for the Rhone River. The NO 3 increase could be explained by the additional inputs from human activities which often dominate the natural sources, and it has been postulated that between 1970 and 1990, humans increased the global delivery of dissolved inorganic N to the oceans by a factor of three , as a consequence of increased consumption of nitrogen fertilizers from 1960 to 1990 by a factor three (from 200 kt.y −1 to more then 600 kt.y −1 ) for the watersheds of the main rivers (cf. Fig. 3). In addition high application rates of manure are found in the northern Mediterranean drainage basins, such as those of the Ebro and Rhone rivers (cf Fig. 4). However, there is an inconsistency between model outputs and observations for recent years, where Ludwig et al. (2009)  discharge at around 40 kt.y −1 at the end of 1990s, while the LPJmL-Med simulates more then 100 kt.y −1 for the same period.

Spatio-temporal variability of nutrients
The rapid increase in NO 3 discharge in the model is clearly associated with the high NO 3 concentration in fertilizer and manure over the watershed of the Ebro, as a consequence of the impressive growth of agrarian production in Spain, multiplied by 3.33 between 1900 and 2008 (Molina et al., 2016).
For the Nile river (Fig. 9d) where nitrate data are missing, LPJmL-Med simulates an increasing discharge of NO 3 from 625 1970, certainly due to the huge amount of NO 3 in the inputs data (cf. section 4) due to the demographic explosion in Egypt and to changes in agricultural practices. In addition the Nile River has one of the highest average nitrate concentrations (Ludwig et al., 2003), although this value has only been derived from a few published values, and it is not clear whether this value is really representative for this river (for more detail see Ludwig et al., 2003).

630
Phosphate (in kt.y −1 of P-PO 4 noted here PO 4 ) output trends are more heterogeneous (Fig. 10). There is a strong increase in the PO 4 discharge simulated for the four main Mediterranean rivers at the beginning of the 1960s. PO 4 discharges reached their maximum value between 1975 and 1980 for the three rivers (Rhone, Po, and Ebro), and relatively later for the Nile river.
This evolution stops after the beginning of 1980s, and the value started to decrease continuously, but the pattern of change is not completely in phase between the main rivers, since the beginning of the phosphate decline differs between the three. It 635 started earlier in the Rhone and Ebro (about 1980), followed by the Po (around 1985), and finally by the Nile (after 1990). This difference reflects the time lag between the different countries in the regulation of phosphorus pollution through the banning of phosphorus detergents and the systematic implantation of wastewater treatment plants. At the end of the 1990s the phosphate discharges again reach the values they had at the beginning of the 1970s'. This indicates that the upgrading of wastewater treatment has been successful. Phosphorous load from industry has also been reduced due to the use of cleaner technology.

640
The comparison with the observed PO 4 from Ludwig et al. (2003Ludwig et al. ( , 2009, and from MOOSE data-set (Raimbault and Lagadec 2012) shows a good agreement with simulated PO 4 , especially for the Rhone River (Fig. 10b) discharge than that from the Ebro River. In the case of the Pô River, it is likely that the high simulated values are associated with the high level of PO4 in fertilizers (cf. Fig. 3d), and in wastewater (cf. Fig.5d). In contrast, PO 4 concentrations are very low in fertilizers and wastewater input for the Ebro River (cf. Fig.3d and Fig. 5d).
in Fig.10e does not show a convex shape like in the three other rivers. A very slow decline of PO 4 is simulated in the Nile River after 1990 as compared to the Europeans countries, thereby indicating that the progression to better water quality has been more rapid in western Europe than in southern countries. Figure 11 shows the NO 3 and PO 4 inputs to the eastern basin (EMed) and western basin (WMed) through river discharge and 655 water runoff transported across the whole basin in the LPJmL-Med grid (cf. Fig. 2). It can first be seen that both NO 3 and PO 4 in water discharges are much higher in the EMed than in the WMed. The differences between WMed and EMed are clearly associated with the particularly high nutrient concentrations in the Nile and Po rivers, which have a considerable impact on the budgets of the entire eastern Mediterranean. Comparison with the in-situ data from Ludwig et al. (2009) suggests a good agreement with the LPJmL-Med outputs for the two basins, the difference between the EMed and the WMed being quite well 660 reproduced by the model with a higher discharge in the EMed (Fig. 11). The huge increase in PO 4 discharge at the beginning of 1960s (from 20 to 80 kt.y −1 in the WMed and from 55 to 160 kt.y −1 in the EMed) is particularly well simulated by the model, as well as the PO 4 drop at the end of the 1980s (Fig. 11b). This huge decrease reflects the adoption of new legislation by the surrounding countries, such as the prohibition of phosphorus detergents, and the improvement of wastewater treatment plants, as well as other features (see below). By contrast, NO 3 discharge shows a constant increase over the whole period (especially 665 in the WMed) which is well reproduced by the model (Fig. 11a), but as already analyzed for the rivers Rhone, Po and Ebro, the model simulates higher NO 3 concentrations during the recent period, particularly for the Ebro River (Fig. 11a, Fig. 9).

Comparison of P and N trends in the Mediterranean Sea
The patterns of change in the river fluxes of N and P exhibit opposite trends in the Mediterranean Sea after the 1980s. Both enhanced by anthropogenic activities in the drainage basins at the beginning of the 1960s, especially P with a dramatic increase in the 1960s and 1970s. However, the efforts undertaken to mitigate point source pollution at the beginning of the 1980s had 670 an immediate impact on the PO 4 loads from rivers, in particular European rivers, where the mean annual fluxes decreased from the second half of the 1980's. This pattern of change in P is mainly due to the improvement of the water quality, with the introduction of tertiary treatment (i.e. with phosphorous removal) and to the banning and abandon of phosphate detergents at the beginning of 1980 especially in Europe. Crouzet et al. (1999) have estimated that between 50 and 75% of the dissolved phosphorous is derived from point sources (i.e. urban wastewater), which confirms the stronger dependence of phosphorous 675 loads on point source pollution, while the other sources (i.e. agriculture) generally account for 20 to 40% (Ludwig et al., 2009).
These trends are in good agreement with the the changes in PO 4 published in Ludwig et al. (2003Ludwig et al. ( , 2009). On the other hand, the pattern is different for NO 3 , which followed a more or less a continuous increase from the beginning of the 1960s until the present in all three rivers. N input is usually dominated by diffuse sources such as fertilizers from agriculture (Ludwig et al., 2003). Crouzet et al. (1999) have estimated that about 45-90% of the nitrogen load in inland waters is related to agriculture in 680 Europe. Furthermore, wastewater and deposition of nitrogen oxides (NOx) in land and surface water may also have contributed to increased N during the recent years (Ludwig et al., 2009).
Finally, a quantitative comparison between in-situ data and model outputs has also been undertaken using a Taylor diagram (Taylor, 2001) that shows the standard deviation, the correlation between the data and the model and the root mean square error for water discharge (Fig. 12a), NO 3 (Fig. 12b) and PO 4 (Fig. 12c). For water discharge, correlations for the three rivers (Rhone,

685
Po, Ebro) are higher than 0.94, with relative standard deviations ranging between 0.7 and 0.8, showing that the LPJmL-Med can predict relatively well the water discharge for the main rivers in the Mediterranean Sea for the historical period (between 1920 and 1980).
For nutrients, the correlation with data is less favorable, and certainly some aspects in the simulation still need to be improved, regarding both the representation of processes and input data (see below). The poor correlation could also be partly 690 ascribed to the limited number of available long time series of in-situ data, which is well below the number of data for the water discharge. The only long-term data presented in Ludwig et al. (2003Ludwig et al. ( , 2009) are produced by the combination of in-situ data and the NEWS-DIN model, which makes more complicated the development of modelling approaches and their evaluation.
In the case of the Rhone river, the MOOSE observatory system provides regular monitoring of nutrient fluxes in the Rhone.
However, in-situ data from the MOOSE program show lower concentrations during the peak of PO 4 in the 1980s as compared 695 to data from Ludwig et al. (2003) (i.e ∼ 9 kt.y −1 in Ludwig et al. (2003) and ∼ 6 kt.y −1 from MOOSE program). The monitoring programs are generally based on monthly/yearly sampling, and therefore not appropriate to quantify the short and violent flash-floods which are typical for the hydrological regimes in the Mediterranean climate (Estrela et al., 2000). If these events are not monitored by suitable sampling strategies, it is not possible to assess the average fluxes, and the ratio of peak discharge to mean annual discharge is still very high. This is particularly true for the long period for the small Mediterranean rivers 700 (Estrela et al., 2000). Furthermore, the Mediterranean river basins are highly reactive to local climatic features, they may be more vulnerable in the context of climate change and alterations in the frequency of extreme climatic events (floods, droughts) can have a severe impact on the river fluxes.
Moreover, the module to simulate nutrient discharge was implemented for the first time in LPJmL-Med model, and it was tested for the first time in this work. Several mismatches between model and observations could be associated with our approach 705 as described above (see Sect. 2), but also to the input data (of fertilizer, wastewater, manure, and land-use patterns used in this study) described in section 2 and 3. In particular, the IFA input time-series considers only national data. Gridded data sets exist, but they are available for specific periods only. Their combination with yearly gridded crop areas and country-level and (when possible) sub-national fertilizer data could enable the establishment of fertilizer input data that better reproduce the spatio-temporal dynamics which could allow a better representation of the spatial differences in fertilizers inputs. The same 710 holds true for manure inputs. Nevertheless, all these comparisons (Fig.9, Fig.10, Fig.11, and Fig.12) led us to conclude that LPJmL-Med was able to reproduce most of the major nutrients (i.e. NO 3 and PO 4 ) features and this provided a basis for realistic values for water and nutrient discharge of the major river into the Mediterranean Sea, as well as at regional scale for the two basins, with a significant correlation and a relatively high R-square for the two basins (ranging between 0.7 and 0.8 for NO 3 , and between 0.6 and 0.7 for PO 4 , Fig. 12). These results suggest that this approach is appropriate for generating a 715 simulation that is sufficiently realistic on decadal timescales, and could be used to investigate the effects of the variations of river discharge in the Mediterranean Sea on marine ecosystems. However, some aspects in the simulation still need to be improved especially the inconsistency between model and in-situ data of NO 3 for recent years: where data remained approximately constant, the model simulates à rapid increase of NO 3 for almost all the rivers. An improved representation of some of the involved processes in future work could improve the simulation of NO 3 in river flow (note that very few in-situ data of NH4 720 are available, that hardly allows us to evaluate the model performance for this state variable). One of the main limitations of the present study is the number of available in-situ data covering the whole Mediterranean basin for input forcing data (i.e. fertilizer, manure, wastewater release data) and for nutrient concentrations in runoff and in river discharge. Therefore, it will be necessary to incorporate new in-situ data (as previously explained), or to develop new statistical approaches to validate our modeling approach in future work.

Conclusions
This study proposes the first basin-wide simulation at 1/12°of water discharge and nutrient release (N and P) into the Mediterranean Sea through the implementation of the biogeochemical land-sea nutrient transfer processes within the agro-ecosystem model LPJmL-Med. For this purpose, the representation of the nutrient transfer from land to sea has been introduced into LPJmL-Med by considering the following processes: remineralization, denitrification, adsorption, nitrification, and phyto-730 plankton dynamics, and a compilation of a new input data set of fertilizer, manure and wastewater nutrients content has been added to the model.
The model successfully simulates the interannual variability of water discharge for the main rivers in the Mediterranean Sea especially the rivers Po, Rhone and Ebro, where we find a very high correlation (the R-square values are higher than 0.94 for the three rivers). The Ebro river show a net decrease of water discharge between 1960 and 1990 (from 28 to 9 km 3 .y −1 ) due 735 to the higher frequency of dry periods observed in the Ebro River basin during the same period. However, the general trend for the Po and the Rhone rivers remains about constant for the same period; this distinctive behaviour pattern may be explained by the non-Mediterranean climate in the upper northern parts of their basins, which results in a high precipitation load. The succession of dry and wet periods is well simulated by the model relative to in-situ data for the three rivers, indicating that the seasonal variability in water discharge is well simulated by the LPJmL-Med. However, the amplitude of the seasonal cycle is 740 greater then observations because we do not represent explicitly the role of dams in the regulation of water flows (i.e., release and retention of water). It is very hard to evaluate and simulate the discharge from the Nile, which loses most of its water through infiltration in swamps, river evaporation and anthropogenic water use. and the construction of the Aswan High Dam in 1965 had a major impact on the water discharge.
Results show a good consistency between the simulated nutrient concentrations (NO 3 and PO 4 ) and available in-situ data.  Code availability. The original code of the LPJmL model is publicly available through PIK's gitlab server at https://gitlab.pik-potsdam.de/lpjml/LPJmL.
The source code of the adjusted LPJmL-Med version developed in this study and described here are available through osupytheas's gitlab server at https://gitlab.osupytheas.fr/mayache/lpjml-med_version1 The output data from the model simulations described in this study can be obtained from the corresponding authoron reasonable request, and should be referenced as Ayache et al.
Data availability. The data associated with the paper are available from the corresponding author upon request. All the data used in this study were published by their authors as cited in the paper. Here we present the model result against the in situ data already published in the literature.    Discharge database RivDIS (Vörösmarty et al., 1996) at the mouths of the a) Rhone, b) Po and c) Ebro rivers.  (Vörösmarty et al., 1996(Vörösmarty et al., ) (averaged over 1920(Vörösmarty et al., -1980.   outputs for the Rhone River (in red), the Po River (in blue) and the Ebro River (in green), western basin (in brown), and eastern basin (in purple) . a) water discharge (data from Vörösmarty et al. 1996Vörösmarty et al. betwwen 1920Vörösmarty et al. and 1985, b) NO3 concentrations (data from Ludwig et al. 2009Ludwig et al. , between 1960Ludwig et al. and 2000, c) PO4 concentrations (data from Ludwig et al. 2009Ludwig et al. , between 1960Ludwig et al. and 2000. RATIO C:P,phyto C:P ratio of phytoplankton and water decomposers gC/gP 250 Elser et al. (2000) sed Fraction of residues trapped daily into sediments d −1 0.024 Billen et al. (1994) VMAX res,water Maximum daily remineralization rate in the water d −1 0.06 Billen et al. (1994)