Hydrological and ecological controls on dissolved carbon concentrations in groundwater and carbon export to surface waters in a temperate pine forest watershed

1 Laboratoire Environnements et Paléoenvironnements Océaniques et Continentaux (EPOC), CNRS, Université de Bordeaux, Allée Geoffroy Saint-Hilaire, 33615 Pessac Cedex France. 2 INRA, UMR 1391 Interactions Sol-Plante-Atmosphère (ISPA), 33140 Villenave d’Ornon, France. 10 3 Departamento de Geoquímica, Universidade Federal Fluminense, Outeiro São João Batista s/n, 24020015, Niterói, RJ, Brazil. *Now also at Laboratoire d’Océanographie et du Climat, Expérimentations et Approches Numériques (LOCEAN), Centre IRD France-Nord, 32, Avenue Henri Varagnat, F-93143 Bondy, France. 15


Introduction
Since the beginning of the Industrial Era, human activities have greatly modified the exchanges of carbon between the atmosphere and continents, as well as those occurring along the aquatic continuum that connects the land and the coastal ocean (Ciais et al., 2013;Regnier et al., 2013).Globally, the land (i.e., vegetation and soil) is a major reservoir of carbon that acts as a net carbon sink (Ciais et al., 2013).However, how and where this carbon is stored in land or exported to inland waters is still a matter of uncertainties (i.e., "the missing anthropogenic CO 2 sink") (Cole et al., 2007).The "missing sink" is the difference between anthropogenic inputs of CO 2 in the atmosphere, the measured increase of CO 2 in the atmosphere and fairly well constrained estimates of the net uptake of CO 2 by the world's ocean.We do not know to what extent inland aquatic systems matter to "missing sink" (Cole et al., 2007).In addition, the amount of C exported from land to streams and rivers is significant compared to the net terrestrial sink and other anthropogenic fluxes (Cole et al., 2007;Ciais et al., 2013;Regnier et al., 2013).The subsidies of carbon from land are also the major source of CO 2 degassing in streams and rivers (Cole et al., 2007;Hotchkiss et al., 2015) and of organic and inorganic C export to the coastal zone (Meybeck, 1982).Hence, it appears crucial to better understand the mechanisms that control storage and export of carbon in land.
Terrestrial vegetation takes carbon from the atmosphere through photosynthesis (GPP, Gross Primary Production), part of which is used by plants as a source of energy and then directly released by autotrophic respiration, while another part is assimilated by vegetation to produce biomass (NPP, Net Primary Production) (Porporato et al., 2003).The C fixed by NPP enters into the soil as organic carbon through mortality and litterfall, root detritus and root exudates (Davidson and Janssens, 2006).The soil carbon is utilized by heterotrophic organisms which respiration is called heterotrophic respiration (Raich and Schlesinger, 1992).Thereafter, soils can lose carbon from export through hydrological processes (Dawson and Smith, 2007).Hydrological processes consist in surface erosion through surface runoff, channel erosion of streams and rivers as well as lateral drainage (Dawson and Smith, 2007).Hydrological C export is generally controlled by basin slope, rainfall intensity, river flow intensity and lithology.Furthermore, soils can lose carbon in groundwaters through vertical leaching.
Thus, drainage of groundwater enriched in dissolved carbon acts as an important source of CO 2 degassing (Venkiteswaran et al., 2014).However, dynamics on how soil carbon interacts with underlying groundwater is poorly documented.
In this study, we focused on dissolved carbon (inorganic [DIC] and organic [DOC]) dynamics in groundwater and associated streams draining a temperate pine forest watershed.We investigated what controls dissolved carbon temporal storage in groundwater and carbon export to superficial waters, including hydrological factors (precipitation, evapotranspiration, water table, drainage) and ecological factors (net ecosystem exchange, gross primary production, ecosystem respiration).We selected as main study site a forest plot, located in a small temperate watershed which offers the convenience of a relatively homogeneous lithology (sand) and vegetation (pine forest), as well as a simple hydrological functioning (outputs mainly as groundwater drainage; no surface runoff).We aimed to describe the factors controlling the leaching of carbon from the soil to groundwater as well as the export of carbon from groundwater to streams.

Study site
The Leyre watershed (2,100 km²) is located in the southwestern part of France near Bordeaux and lies in the "Landes de Gascogne" area (Fig. 1).The landscape is a very flat coastal plain in almost all its surface with a mean slope lower than 1.25 ‰ (generally NW-SE), but with local gentle slopes (notably near some streams) (Jolivet et al., 2007).The mean altitude is lower than 50 m (Jolivet et al., 2007).The lithology is relatively homogeneous and constituted of sandy permeable surface layers dating from the Plio-quaternary period (Legigan, 1979).The soil is composed of sandy permeable podzols characterized by a low pH (4), low nutrient availability, and high organic carbon content that can reach 55 g per kg of soil (Augusto et al., 2010).The region was a vast wetland until the XIXth century, when a wide forest of maritime pine (Pinus pinaster) was sown following landscape drainage from 1850.Nowadays, the catchment is mainly occupied by pine forest (about 80%), with a modest proportion of croplands (about 15%) (Fig. 1).The typical rotation period of pine forest is 40 years, ending in clear-cutting, tilling and re-planting (Kowalski et al., 2003).The climate is oceanic with mean annual air temperature of 13°C and mean annual precipitation of 930 mm (Moreaux et al., 2011).Moreover, the average annual evapotranspiration is in the range of 234-570 and 63-800 mm, respectively for maritime pine and cropland (Govind et al., 2012).Hence, owing to the low slope (i.e., < 1.25 ‰) and the high permeability of the soil (i.e., hydraulic conductivity is about 10 -4 m s -1 , Corbier et al., 2010), surface runoff cannot take place in the Leyre watershed, and thus the excess of rainfall percolates into the soil and supports the enrichment of carbon in groundwater.In addition, the sandy permeable surface layers contain a free and continuous water table strongly interconnected with the superficial river network.This interconnection is facilitated by a dense network of drainage ditches, initiated in the XIX th century, and currently maintained by forest managers in order to increase tree growth rate.Consequently, hydrology in the Leyre watershed occurs almost exclusively through drainage of groundwater.Furthermore, the seasonal changes in groundwater table can be important, with a water table close to the surface during winters and levelling down to 2.0 m depth below the surface during most summers.
At the Bilos site, the impermeable layer is located at approximatively 10 m below the soil surface.Finally, we adapted the Strahler definition of first order stream by including streams and ditches either having no tributaries and/or being seasonally dry.
All the abbreviations used in this paper are listed in table 1. Carbon fluxes are in italics.

Sampling strategy
Within the Leyre watershed, we selected 3 piezometers that were located in different forest plots and 6 first order streams whose watersheds were dominated largely with forest (~90 %) (Fig. 1).
From 2001, the Bilos site has been equipped with continuous measurements for groundwater table depth, atmospheric fluxes of energy, water vapor and CO 2 and ancillary variables.In addition to these continuous measurements, we monitored the partial pressure of CO 2 (pCO 2 ), total alkalinity (TA) and dissolved organic carbon (DOC) in groundwater and in first order streams (Deirmendjian, 2016).The Bilos site groundwater was sampled with a frequency of approximately once a month, on 15 occasions between Feb-2014 and Jul-2015.In addition, the two other piezometers were sampled respectively on 11 (Aug. 2014-Jul. 2015) and 6 occasions (Jan. 2015-Jul. 2015).The first order streams were sampled on 16 occasions between Jan. 2014 and Jul.2015.Concerning river discharge, our study took benefit from one calibrated gauging station of the water quality agency (with a daily temporal resolution), located on the Leyre River (4 th order stream) (Fig. 1).

Continuous measurements at the Bilos site
Precipitation was measured using automatic rain gauges with a 30 minutes integration: one Young EML SBS 500 (EML, North Shields, UK) was located in a small clear-cut at 3 m above ground from 01/01/2014 to 10/05/2014 and one electronic gravimetric heated precipitation gauge TRwS (MPS system; Bratislava, Slovakia) was located at the top of the canopy on a 6 m tower, from 01/07/2014 to 31/12/2015.Hence, between 11/05/2014 and 31/06/2014, there was no precipitation measurement available at the Bilos site.Thus, during this period, we used data from Meteo France © station at Belin-Béliet (about 30 km from the Bilos site).Precipitation measurements were also checked weekly in the field with manual reports.
The net CO 2 and latent heat fluxes were measured using the eddy covariance technique.The eddy covariance technique allows determine the turbulent-scale covariance between vertical wind velocity and the scalar concentration of sensible heat, CO 2 or H 2 O, measured near the ecosystem/atmosphere interface which is an atmospheric flux between the ecosystem and atmosphere.The atmospheric exchange originates from atmospheric eddies (turbulence) caused by buoyancy and shear of upward and downward moving air that transport gases such as CO 2 and H 2 O.
Here, wind velocity, temperature and CO 2 /water vapor fluctuations were measured with, respectively, a sonic anemometer (model R3, Gill instruments Lymington, UK) and an open path dual CO 2 /H 2 O infrared gas analyzer (model Li7500, LiCor, Lincoln, USA) at the top of a 9.6 m tower (01/01/2014 to 10/05/2014) and with another sonic anemometer (model HS50, Gill instruments) and a close path enclosed dual CO 2 /H 2 O infrared gas analyzer (model Li7200, LiCor ©) at the top of a 15 m tower (09/07/14 to 31/12/2015) .Consequently, there were no eddy covariance measurements available between 11/05/2014 and 08/07/2014 and thus between these two dates the latent heat fluxes were determined following the procedure of Thornthwaite (1948).
In this paper raw data were processed following the ICOS methodology (Aubinet et al. (1999)).The post-processing software EddyPro (www.licor.com)was used to treat raw data and compute average fluxes (30 min period) by applying the following steps: (1) spike removal in anemometer or gas analyzer data by statistical analysis, (2) coordinating rotation to Biogeosciences Discuss., doi:10.5194/bg-2017-90,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.align coordinate system with the stream lines of the 30 min averages, (3) linear de-trending of sonic temperature, H 2 O and CO 2 channels (4) determining time lag values for H 2 O and CO 2 channels using a cross-correlation procedure, (5) computing mean values, turbulent fluxes and characteristic parameters, (6) high frequency corrections via transfer functions (Moore, 1986) and ( 7) performing a Webb Pearman Leuning correction to account for the effects of temperature and water vapor on measured fluctuations in CO 2 and H 2 O (Webb et al., 1980).Thereafter, CO 2 and H 2 O fluxes were filtered in order to remove points corresponding to technical problems, meteorological conditions not satisfying eddy correlation theory or data out of realistic bounds.Different statistical tests were applied for this filtering: (1) stationarity and turbulent conditions were tested with the steady state test and the turbulence characteristic test recommended by Foken and Wichura (1996) and Kaimal and Finnigan (1994).
Based on several tests, only values of CO 2 and H 2 O fluxes that pass all the filters mentioned above were retained.Then, missing values of CO 2 and H 2 O fluxes were gap-filled and partitioned into GPP and R eco with the R package Reddyproc (version 0.8-2) that implements the procedure of Reichstein et al (2005).NEE was then partitioned into GPP and R eco by applying the following steps: (1) during nighttime GPP=0 so NEE = R eco (2) statistical regression between R eco and night air temperature and meteorological conditions is adjusted with a Arrhenius type equation (Lloyd and Taylor, 1994) : Where, T soil is the soil temperature measured at 10 cm.T ref , R 10 , E 0 and T 0 are respectively a temperature of 283.15 K, the ecosystem respiration for a reference soil temperature of 10 °C, the activation energy and a calibrated temperature (227.13K).
(3) day-time R eco is obtained by extrapolating night-time fluxes using the temperature response (4) GPP is calculated as the difference between daytime NEE and R eco , additional checks are performed to avoid unrealistic values of GPP.
GPP is positive or zero.R eco is positive.NEE = R eco -GPP.Hence, positive NEE indicates an upward flux whereas a negative NEE indicates a downward flux.We also considered that NEE = F c .
The groundwater table was measured using high performance level pressure sensors (PDCR/PTX 1830, Druck and CS451451, Campbell Scientific) in one piezometer located amid the Bilos site.The pressure measurements were fully compensated for temperature and air pressure fluctuations.The measurements were obtained at 60 second intervals and integrated on 30 min period.They were checked with manual probe weekly.Since there were no measurement available between 30/04/2014 and 23/06/2014, values of the groundwater table depth was interpolated between these two dates.
Further away, we used the parameter H m that is the mean groundwater table between two sampling dates.

Discrete sampling
We measured partial pressure of CO 2 (pCO 2 ) directly in the field and total alkalinity (TA) and dissolved organic carbon (DOC) back in the laboratory.
Thus, partial pressure of CO 2 in groundwater and in first order streams was measured directly using an equilibrator.This equilibrator was connected to an Infra-Red Gas Analyzer (LI-COR®, LI-820), which was calibrated on the 0-90,000 ppmv range, following the procedure of Deirmendjian (2016).We took the precaution to renew the water in the piezometers by pumping of about 300 L with a submersible pump before sampling.
TA was analyzed on filtered samples by automated electro-titration on 50 mL filtered samples with 0.1N HCl as titrant.
Equivalence point was determined with a Gran method from pH between 4 and 3. Precision based on replicate analyses was better than ± 5 µM.For samples with a very low pH (<4.5), we bubbled the water with atmospheric air in order to degas CO 2 .Consequently, the initial pH increased above the value of 5, and TA titration could be performed (Abril et al., 2015).
We calculated dissolved inorganic carbon (DIC) from pCO 2 , TA, and temperature measurements using carbonic acid dissociation constants of Millero (1979) and the CO 2 solubility from Weiss (1974) as implemented in the CO 2 SYS program.
Contrary to pCO 2 calculation from pH and TA (Abril et al., 2015), DIC calculation from measured pCO 2 and TA was weakly affected by the presence of organic alkalinity, because 80±20 % of DIC in our samples was dissolved CO 2 .
DOC samples were obtained after filtration, in the field through pre-combusted GF/F filters (porosity of 0.7 µm); DOC filtrates were stored in pre-combusted Pyrex vials (25 mL) and acidified with 50 µL of HCl 37 % to reach pH 2 and kept at 4 °C in the laboratory before analysis.DOC concentrations were measured with a SHIMADZU TOC 500 analyzer (in TOC-IC mode), which is based on thermal oxidation after a DIC removal step (Sharp, 1993).The precision (repeatability) was better than 0.1 mg L -1 .

Water balance at the Bilos site
By definition, the water balance equation is as follows: Where, P, D, ETR, GWS and ΔS are respectively, precipitation, drainage, evapotranspiration, groundwater storage and change of soil water content in the unsaturated zone, all expressed in mm d -1 .These 5 parameters were determined respectively as follows: (1) P is the cumulative precipitation over a period t as measured at the site; (2) D is the drainage at the Bilos site, estimated as the mean Leyre River flow over a period t divided by the catchment size at the gauging station; (3) ETR is the cumulative evapotranspiration obtained from eddy covariance measurements of latent heat flux over a period t; Biogeosciences Discuss., doi:10.5194/bg-2017-90,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.
(4) GWS is the groundwater storage estimated as the net change in water table depth over the period t times the soil effective porosity; (5) S.No reliable measurement of soil water content was available and this term was not measured therefore.
In the above description, drainage (D) is calculated for the Leyre River that is fourth order stream according to Strahler classification.The Bilos site is drained in majority by ditches that are first order stream according to Strahler classification.
Hence, we corrected the drainage of fourth order stream to estimate the drainage value at the Bilos site.This correction factor is about 2.3, as deduced from in situ discharge measurements in the Leyre watershed (Deirmendjian., 2016).
Finally, water balance of the Bilos site, was calculated on a monthly basis over a two years period (2014)(2015) and also between each measurement period.

Groundwater carbon fluxes at the Bilos site
In order to understand the dynamics of carbon in groundwater, we calculated 4 different terms of carbon groundwater fluxes at the Bilos site: storage of DIC and DOC (DIC s and DOC s ) and export of DIC and DOC (DIC ex and DOC ex ) all expressed in mmol m -2 d -1 .
Storage of DIC in groundwater is calculated using the following equation: Where, S f and S i are the final and the initial stock of DIC in groundwater in mmol m -2 .DIC f and DIC i are the final and the initial concentration of DIC in groundwater in mmol m -3 , respectively.V f and V i are the final and the initial volume of groundwater in m 3 m -2 .dt is the period in day between two sampling days.DIC (or DOC) storage can be positive or negative depending if gain or loss of DIC occurred in the groundwater between two sampling days.
The volume of groundwater (V) was calculated as the following manner: Where, 10 and H (H is negative), are respectively the total height of the permeable surface layer and the height of groundwater table, in m.In the Leyre watershed the total porosity equals to 0.4 whereas Φ effective (effective porosity) equals to 0.2 (Augusto et al., 2010;Moreaux, 2012) We calculated DOC s as the same manner as DIC s .
Export of DIC through drainage of groundwater is calculated using the following equation: Where, S DIC is the mean stock of groundwater DIC between two sampling dates in mmol m².S f and S i are the final and the initial stock of DIC in groundwater in mmol m -2 .DICi and DIC f are the initial and the final concentration of DIC in groundwater in mmol m -3 , respectively.V i and V f are the initial and the final volume of groundwater in m 3 m -2 .
We calculated DOC rs as the same manner as DIC rs.

Degassing in first order streams
We determined degassing flux (F Degass ) between two sampling dates, in first order streams as described in Deirmendjian (2016).
Where, ΔCO 2 (t) and ΔCO 2 (t+1) are the differences between the concentrations of CO 2 in Bilos groundwater and in first order streams (mean of the 6 first order streams) at time t and t+1, expressed in mmol m -3 .Q mean is the mean river flow of first order streams in m 3 d -1 between time t and t+1.S is the area of the Leyre watershed in m².

Analysis of data
In this paper we used Pearson's correlation coefficient (r p ) to investigate the strength of a linear correlation between mean carbon concentrations in groundwater (DIC m , DOC m ) and in first order streams (DIC m1 , DOC m1 ) with carbon groundwater fluxes (DIC s , DOC s , DIC ex , DOC ex ), hydrological parameters (P, GWS, ETR, D and H m ) ecological parameters (NEE, GPP, R) and degassing flux in first order streams (F Degass ).
All the Pearson's correlation coefficient, r p are presented in Table 2.

Water mass balance
Over the sampling period ( 2014 During the 2014-2015 period, the Leyre River discharge was on average 17.9 m 3 s -1 including two relatively short flood events (further referred as high flow periods) in Jan. 2014-Apr.2014 (maximum flow of 120 m 3 s -1 ) and in Feb. 2015-Mar.
Furthermore, there was a good linear relationship between GWS and P (r p = 0.76, p-value < 0.05), and between GWS and ETR (r p = -0.73,p-value < 0.05) (Tab.2).High river discharge periods were also associated with the highest water table and the highest D (Fig. 2a-b).As a consequence, H m and D were positively correlated (r p = 0.86, p-value < 0.001) (Tab.2).
Overall, ETR was higher in spring and summer (maximum value of 5.33 mm d -1 in Apr.2014) than in autumn and winter (minimum value of 0.33 mm d -1 in Dec. 2014) (Fig. 2b).On the contrary, P was higher in autumn and winter (maximum value of 7.99 mm d -1 in Jan. 2014) than in spring and summer (minimum value of 0.18 mm d -1 in Sep.2014) (Fig 2b).
Finally, water balance at the Bilos site as calculated as P, on the one hand and as the sum of ETR, GWS and D on the other hand closely followed the 1:1 Line (Fig. 3).The water mass balance estimated with different techniques and independent devices was thus fairly consistent.This reveals that our approach was well adapted to establish the water mass balance of our forest plot and thus the dissolved carbon fluxes

Carbon fluxes
In groundwater and first order streams, TA originated from weathering of silicate minerals with vegetation-derived CO 2 (Polsenaere and Abril 2012; Deirmendjian 2016).In addition, the proportion of TA in the DIC pool was respectively 5±5 % and 30±15 % for groundwater and first order streams, the large majority of the DIC was thus composed of dissolved CO 2 resulting from microbial and plant root respiration in the soil.
In high water table and high D periods (Fig. 2a; 4).Furthermore, during these high flow periods, we never observed such high DOC concentrations in first order streams as in the groundwater (Fig. 4c).Thus, DOC m1 was not correlated with DOC m (r p = 0.39) and DOC ex (r p = 0.30) (Tab.2).On the contrary, highest values of DIC in Bilos groundwater (5,400 and 5,100 mmol m -  2a, 4).Moreover, at the same time intervals, we observed a DIC s increase in the groundwater (33 and 52 mmol m -2 d -1 ) roughly equivalent to a DOC s decrease (-31 and -72 mmol m -2 d -1 ) (Fig. 5b).We observed a fast second increase (Aug.2014-Sep.2014) of DIC in Bilos groundwater (2,700 to 5,400 mmol m -3 ) but not related with any decrease of groundwater DOC (Fig. 4b-c).Same temporal trend was observed for piezometer 2 (Fig. 4b-c).
DIC concentration in Bilos groundwater decreased from 5,400 mmol m -3 (Sep.2014) to 1,700 (Mar.2015) mmol m -3 in parallel with a rise in the water table due to high P (Fig. 2b; 4b).Same temporal trend was observed for piezometer 2 (Fig. 2b, 4b).Concomitantly, a fast increase in DOC concentration from 575 to 3,670 mmol m -3 occurred in Bilos groundwater between Jan. 2015 and Mar.2015 (Fig. 4c).Furthermore, the same trend (fast increase of DOC) was observed for piezometer 3 but not for piezometer 2 (Fig. 4c).Then, from Mar. 2015 to Jul. 2015 a large decrease of DOC concentration from 3,670 to 320 mmol m -3 in Bilos groundwater was observed in parallel with a small increase of DIC from 1,700 to 2,400 mmol m -3 , and a drop in the water table (Fig. 2a, 4b-c).During this period, the same trend was observed in piezometer 3 but not in piezometer 2 (Fig. 2a, 4b-c).
DIC m and DOC m were negatively correlated in Bilos groundwater (r p = -0.64,p-value < 0.05) (Tab.2; Fig. 4a).Moreover, DIC m in Bilos groundwater was negatively correlated with D (r p = -0.70,p-value < 0.05) and H m (r p = -0.82,p-value < 0.001) whereas DOC m in Bilos groundwater was positively correlated with D (r p = 0.96, p-value < 0.001) and H m (r p = 0.88, pvalue < 0.05) (Tab.2, Fig. 4a).In addition, DIC and DOC were respectively 2,650±1,270mmol m -3 and 1,250±1,170 mmol m -3 .In the same time, H m controlled both the export of both carbon forms, being positively correlated with DIC ex (r p = 0.83, p-value < 0.001) and DOC ex (r p = 0.77, p-value < 0.05).Although seasonal differences occurred between both carbon forms throughout the sampling period, the mean, time-integrated value of carbon export was 0.86 mmol m - Overall, throughout sampling period, storage of carbon in Bilos groundwater was highly variable, depending on the intensity of increase/decrease of carbon concentrations in groundwater, with mean value of 0.85 mmol m -2 d -1 and -9.6 mmol m -2 d -1 , for DIC s and DOC s respectively (Fig. 5b).This means that throughout sampling period, the groundwater gained some DIC but lost some DOC.Moreover, DIC s was correlated with none of the studied parameters whereas DOC s was correlated with P (r p = 0.63, p-value < 0.05) and GWS (r p = 0.62, p-value < 0.05) (Tab.2).
Mean GPP, mean R eco and mean NEE were respectively 420 mmol m -2 d -1 , 310 mmol m -2 d -1 and -110 mmol m -2 d -1 throughout the sampling period (here we excluded 16/05/14-07/07/14 period, because there were no data available), equivalent to 1,845; 1,355 and 495 g C m -2 yr -1 .This figure is close from Moreaux et al (2011) estimates of 1,720; 1,480 and 340 g C m -2 yr -1 respectively, as measured at a younger forest stage in the same forest plot.NEE was positive (R eco > GPP) in Oct. Nov and Dec. 2014 and negative (R eco < GPP) all over the rest of the sampling period (Fig. 5a).GPP increased from Our dataset -obtained during a 15 months long monitoring-enabled to understand well how rainfall is partitioned between evapotranspiration, drainage and groundwater storage in the "Landes de Gascogne" area, as well as what controls the water table level between precipitation, drainage and evapotranspiration.Due to the high permeability of the sandy podzolic soil, when precipitations are high, water infiltration in the soil is faster than the capture by vegetation.Consequently, groundwater is filled directly by rain water, which raises the water table and thus increases the GWS.For that reason, GWS tightly increases with P (Fig. 2b; Tab. 2).On the contrary, the transpiration flow through plants and the ETR are maximum in spring and summer when the precipitations are minimal.For that reason, GWS decreases with increasing ETR (Fig. 2b; Tab. 2), revealing that water uptake by the pine trees exerts a direct influence on water table depth.
On the one hand, this was consistent with some authors who found that water table was significantly elevated after harvesting pine forest due to reduced evapotranspiration (Bosch and Hewlett, 1982;Sun et al., 2000;Xu et al., 2002;Guillot et al., 2010), and Guillot et al (2010) who estimated that groundwater contribution to the ETR was above 50 %.On the other hand, rising water table may saturate the plant rooting zone, and the putative anoxia may place the vegetation under stress because of lacking oxygen required for aerobic respiration (Naumburg et al., 2005).However, the network of drainage ditches created by foresters evacuates very rapidly the water in excess of field capacity when the groundwater level rises above 0.5 m depth.Since most pine roots are located above this level (Bakker et al., 2006), the pine trees do not exhibit any transpiration reduction when the groundwater is high, as reported by e.g.Loustau et al (1990).In other environment (Mediterranean), spring growth flush (as well as groundwater uptake) of cork oak was also initiated when groundwater was near the soil surface (Costa et al., 2003).In semi-arid oak savanna in California Miller et al. (2010) also showed that woody vegetation uses a significant amount of groundwater as soil moisture reserves are depleted.Moreover, the precipitation preceding the growing season, can be important to simulate physiological activity of the trees during the growing season (Miller et al., 2010;Naumburg et al., 2005).In the studied watershed, climate is oceanic and the precipitation preceding the growing season is always intense (mean precipitation during Nov. Dec.Jan and Feb 2000-2015 was 320 mm).Hence, groundwater uptake by plants in spring could be stimulated by high precipitation in winter and high water table period that enhance soil moisture reserve.
Throughout the sampling period, H m fluctuations control the intensity of D (Tab.2).This correlation linking D and H m is not unexpected since H m might be interpreted as a proxy for the watershed pressure head driving the drainage D (Tab.2; Fig. 2ab).Conversely, H m is not correlated with P, ETR or GWS.
Consequently, the drainage (D) was not correlated with GWS (Tab.2), and fluctuation of groundwater storage was decoupled from D in this ecosystem.This was due to the long time needed for transferring groundwater discharge to river flow in this system.This was indeed consistent with the flat topography of the watershed (mean slope lower than 1.25 ‰) Biogeosciences Discuss., doi:10.5194/bg-2017-90,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.
and its low overall hydraulic conductivity.It should be also noticed that the human-made network of ditches probably has a substantial role in the drainage regime.Indeed, when the water table is high, it is strongly regulated by drainage ditches which are connected to rivers.Finally, fluctuation of the height of the water table was not driven specially by one of the hydrological parameters, but influenced by the water balance as P -ETR-GWS-D.
We consider water input as precipitation and water output as evapotranspiration plus drainage.Water stock changes (i.e., groundwater storage plus soil water in the unsaturated zone) can be calculated from subtracting output to input.Hence, water stock changes were -1 mm and -48 mm, respectively in 2014 and 2015 (Tab.3).However, we measured groundwater storage of -72 mm and -174 mm, respectively in 2014 and 2015 (Tab.3).Thus, we assume that soil water in the unsaturated zone gains 71 mm and 126 mm; respectively in 2014 and 2015.Larger increase of soil water content in 2015 is consistent with the larger thickness of the unsaturated zone in 2015 (up to 1,800 mm) (Fig. 2a).At the Bilos site, Moreaux et al ( 2011) measured a loss of soil water of -70 mm, but calculated on the 0-80 cm soil layer.This suggests that the unsaturated zone of the soil gains water in deeper horizons and lose water in superficial horizon.

Production and consumption of dissolved carbon in groundwaters
DIC in the groundwater was composed of 5±5 % (range is 0-24 %) of bicarbonate (alkalinity) and 97±3 % (range is 76-100 %) of dissolved CO 2 .Dissolved CO 2 originates from the solubilization of CO 2 that comes from soil respiration, root respiration (autotrophic) and microbial respiration (heterotrophic) (Raich and Schlesinger, 1992); Dissolved CO 2 can also originate from respiration in the saturated soil, i.e. in the groundwater itself using DOC (Craft et al., 2002).Total alkalinity (TA) in the sampled groundwater originates from silicate weathering, in the sandy podzols (Polsenaere and Abril, 2012;Polseneare et al. 2013;Deirmendjian 2016).Indeed, dissolved CO 2 can react with silicates to produce bicarbonates However, in monolithic watersheds draining mainly silicate rocks, TA is typically very low, on average below 125 mmol m -3 according to Meybeck (1987), and 32-135 mmol m -3 in groundwater at the three sampled sites.This is consistent with the very low content in feldspars and allover clay minerals in our sandy podzols (Augusto et al., 2010).Moreover, silicate weathering and soil respiration occur at different time scales.Indeed, soil respiration is a short time scale process (̴ 10 0 -10 2 years) whereas silicate weathering is a long time scale process (̴ 10 4 -10 6 years) (Ciais et al., 2013;Colbourn et al., 2015).DOC in groundwater generally comes from the leaching of soil organic matter.High soil pH or calcite-rich soils containing clay favor DOC stabilization, while low soil pH in combination with sandy soils, as the case here in the Landes de Gascogne, favor DOC destabilization and lixiviation (Paradelo et al., 2015).Moreover, in the Leyre watershed, sandy podzols contain almost no clay minerals (Augusto et al., 2010) and, therefore, this absence of phyllosilicates likely prevents the formation of clay-OM complex from occurring, and thus probably slows DOC stabilization in soil.
During our sampling, the temporal variations of DOC and DIC concentrations in groundwater were opposite (Fig. 4), showing that this two carbon forms originate from very different processes occurring in the saturated and/or unsaturated soil.
During flood peaks (high flow periods), water table at the Bios site was high (as a consequence of high precipitation and low vegetation use during winter, see section 4.1) and was very close to the soil surface: 1.2 cm in Feb. 2014 and 17.2 cm in Mar.2015 (Fig. 2a).Consequently, the groundwater had reached the superficial and organic-rich horizons of the soil (0-20 cm) where soil organic carbon concentration is highest, 18.3 g kg -1 in the Bilos site (see section 2.2).There were non-negligible amounts of organic carbon in deeper soil layers, but content values were much lower (6.4 g kg -1 in the 40-60 cm layer) than in the topsoil layer.In addition, organic matter in deep soil layers is generally well-stabilized (Torn et al., 1997;Rumpel and Kögel-Knabner, 2011) and consequently not prone to produce DOC.In their review of stream composition in temperate forests, Michalzik et al. (2001) 2014) in soil waters of forest top layers).Moreover, we calculated a stock of SOC in the 0-60 cm layer of 9.7 kg m -2 , whereas the stock of groundwater DOC ranges from 0.08 kg m -2 (during high flow period) to 0.01 kg m -2 (during base flow period).Finally, only a tiny part of the soil organic matter content is likely to be leached into groundwater.
Furthermore, as storage of DOC in Bilos groundwater increased with GWS (r p = 0.62) and P (r p = 0.63), groundwater gets enriched in DOC when GWS increases due to high precipitation (Tab.2).However, during high flow periods of 2015, DOC in piezometer 2 and 3 did not reach values as high as the Bilos site.This is a consequence of both local hydrological and topographic heterogeneity.Indeed, the piezometer 2 is located very close to a stream, in a gentle slope part, and thus water table level remains lower (maximum level of groundwater table is 157 cm below the soil surface) than in other piezometers and do not reach organic horizons.The piezometer 3 is located on a small forest plot similar to the Bilos site.However, the groundwater table reaches only 70 cm below the soil surface.In addition, at 70 cm depth, DOC concentrations in forested acidic soil are close to 800 mmol m -3 (10 mg L -1 ) (Camino-Serrano et al., 2014).
Immediately after the maximum flow, DOC concentration in groundwater decreased quickly and is associated in parallel with an increase in DIC (Fig. 4b-c).Furthermore, during this DOC decrease, storage of DIC and DOC in groundwater were almost equivalent but opposite: +90 mmol m -2 d -1 for DIC versus -100 mmol m -2 d -1 for DOC between 12/02/2014 and 16/05/2014 (Fig. 5b).In addition, during high flow periods the mean residence time of groundwater DOC was 80±50 days (Fig. 6).This residence time is long enough to assume that the increase of DIC was due mainly to the respiration of groundwater DOC.The respiration rates of 93 mmol m -2 d -1 (i.e., DOC s -DOC ex between 12/02/14 and 16/05/14) in the Leyre watershed is coherent with respiration rates in streams determined by Battin et al. (2008) (i.e., 1.93  molecular weight compounds (Aravena et al., 2004) that also derived from organic matter recently produced and leached.
Thus, during high flow periods groundwater DIC originates mostly from respiration of young labile groundwater DOC.
During base flow periods, groundwater DOC was very stable (535±80 mmol m -3 ) (Fig. 4c), suggesting that groundwater DOC was more recalcitrant and probably more stabilized and more aged during base flow periods.This was consistent with findings of Schiff et al. (1997) which found that a small temperate basin had wide range in DO 14 C, from old groundwater values at base flow under dry basin conditions to relative modern values during high flow or wetter conditions.In the Leyre watershed we found that recalcitrant groundwater DOC has residence time around 30 years in autumn (Fig. 6).In addition, as we never observed an increase of groundwater DOC concentration during base flow periods (when groundwater table is low) (Fig. 4c), it seems that soil DOC in upper horizons cannot be mobilized by rainwater percolation and that saturation of the soil horizon is necessary; this is also attested by the absence of correlation between P and DOC m (Tab.2).However, in other types of environment (fractured rock aquifer) Shen et al. (2015) found significant correlation between surface precipitation and groundwater DOC.
The second increase of DIC during Sep.2014 was due to another process, not associated with any DOC degradation.This DIC likely originated from the dissolution of gaseous CO 2 that could have accumulated in the unsaturated region of the soil during summer and enhanced heterotrophic and plant roots respiration.During drought period (rainfall was 5 mm in Sep. 2014) the water deficits stress growing vegetation and leads to numerous physiological changes (Naumburg et al., 2005).
The dehydration of plants lowers the rate of photosynthesis (1) directly by closing stomatal pores, hence interfering with uptake of carbon dioxide by leaves, and (2) indirectly by adversely influencing the photosynthetic mechanism (Kozlowski, 2002).The second increase of groundwater DIC occurred exactly during positive NEE and larger ecosystem respiration than Gross Primary Production, which we believed is influenced by drought period (Fig. 5a).Thus, when the forest ecosystem is a source of CO 2 for the atmosphere, it is also a source of CO 2 for the underlying groundwater.This is also suggested by the positive correlation (r p = 0.48) between NEE and DIC m (Tab.2).
Moreover, mature forests are in general characterized by a thick humus layer that could in part isolate soil air from atmospheric air and participate to the process of soil CO 2 accumulation during summer.However it is not the case here because the Bilos site is a young forest, trees were sowed in 2005 and after tillage (i.e., tillage buried the humus layer).Such transfer of CO 2 from soil air to groundwater is typical of events of percolation of rainwater in the unsaturated soil after a long dry period as reported in an Amazonian plot (Johnson et al., 2008).However, at our study site, Sep.2014 was one of the driest months throughout the sampling period (Fig. 2b), which suggests that soil CO 2 could have been transported by simple downward diffusion.Thereafter, in Oct. 2014 when the water table started to increase again, groundwater DIC decrease as a consequence of dilution with rainwater with low DIC content.Moreover during this period of rising water table, draining of groundwater was stronger, which resulted in a faster recycling of the DIC in the groundwater (Fig. 2a-b; 4b).Indeed, the calculated mean residence time of groundwater DIC was 80±60 days during the Oct. 2014-Dec.2014 period, while it was 3,200±3,000 days during the other periods.(Fig. 6).Furthermore, residence time of DIC and DOC (i.e., calculated between two sampling days) ranged from months to several years (Fig. 6).This reveals the intensity of different

Carbon export from groundwater and fate in streams
In the groundwater, the water table (H m ) controls at the same time the intensity of the drainage D, and the concentrations of dissolved carbon, positively for DOC and negatively for DIC (Tab.2; Fig. 4a).In addition, H m and D control the export of DOC and DIC through groundwater.Indeed, when H m is high, D is high and thus DIC ex and DOC ex are both maximal (Tab. 2).Thus, drainage intensity and DOC concentration in groundwater have a cumulative positive effect on DOC export; in contrast, drainage intensity and DIC concentration in groundwater have antagonist effects on DIC export, but, as the drainage effect is stronger, DIC export is still positively correlated with H m and D. As a result, throughout the sampling period, high flow and high H m periods (Jan. 2014-Apr. 2014and Feb. 2015-Mar. 2015) account for 90 % and 50 % of the total DOC and total DIC exports, respectively.
Consequently, groundwater exports the majority of DOC during high flow periods, but about the same quantity of DIC during base flow periods and high flow periods.However, for the whole sampling period, the mean carbon export is almost the same (about 0.85 mmol m -2 d -1 ), both for DIC and DOC (Tab.4), and, as drainage is the only hydrological pathway, the forest ecosystem exports in total 1.70 mmol m -2 d -1 , 50% as DOC and 50% as DIC.Throughout the sampling period, this export flux represents 1.5 % of the mean NEE (Tab.4).Hence, young forest at the Bilos site, loss a low quantity of carbon through vertical soil leaching and groundwater drainage.In contrast, in peatland systems, Billett et al. (2004) showed that the amount of carbon exported in surface waters can potentially exceed NEE.It should be also noticed that another part of the NEE is lost by root exudates, litter fall and fine roots turnover (Davidson and Janssens, 2006).Afterwards, part of the litter fall and fine roots turnover is respired by heterotrophic respiration in the soil (Raich and Schlesinger, 1992) while another part prone to produce SOC that can be more or less stable (Rasse et al., 2005).
In first order streams, DIC and DOC concentrations were generally lower than in groundwaters, but still negatively correlated, with maximum DOC at high flow and maximum DIC at base flow (Tab.2; Fig. 4b-c).
DOC m1 is controlled positively by H m (Tab.2).Indeed, increase in riverine DOC concentrations with discharge and high water table has been reported in the Leyre watershed and in many other catchments elsewhere.However, during high flow periods we never observed DOC concentration in first order streams as high as those in Bilos groundwater (Fig. 4c).This might be due partly because groundwater DOC is quickly respired before reaching the stream (see section 4.2).Also, as groundwater DOC enters the superficial river network through drainage it might be rapidly recycled by photo-oxidation or respiration within the stream, as attested by the absence of correlation between DOC m and DOC m1 .Indeed, groundwater DOC during high flow periods must be relatively labile (Aravena et al., 2004).In addition, podzols contain a huge quantity of humic and fulvic acids (Righi and Wilbert, 1984), that are also very photodegradable (Tranvik, 1996;Schmitt-Kopplin et al., 1998;Suhett et al., 2007).
The comparison of DOC concentration between the three piezometer (Fig. 4c) also suggest some significant heterogeneity in DOC export due to topographic effects and the depth of the water table, the Bilos plot potentially exporting more DOC than the entire Leyre watershed.As a consequence, part of the difference in DOC concentration between groundwaters and streams during the high flow might also be due to spatial mixing of groundwaters from different soil horizons.In contrast, during base flow period, DOC concentrations in groundwaters and in first order streams were very similar (Fig. 4c), which suggests that this DOC was not labile, and not degraded in the superficial river network.
DIC concentration in first order streams increased in parallel with those in groundwater, during base flow (Fig. 4b).Indeed, concentrations of DIC show an inverse relationship with discharge in many catchments, as the result of dilution with rain water and lower contribution of deep CO 2 -enriched groundwater during high flow periods (Billett et al., 2004;Dawson and Smith, 2007).In addition, during both base flow and high flow periods, we never observed DIC concentration in first order streams close to groundwater DIC concentration.Indeed, the quick loss of DIC between groundwater and first order streams is due to the degassing of CO 2 (F Degass ) (Venkiteswaran et al., 2014).This rapid degassing was also attested by the change in the δ 13 C signature of the DIC (Deirmendjian, 2016).Furthermore, the positive correlation between F Degass and both DIC ex and H m reveals that groundwater DIC is the main source of CO 2 degassing in superficial stream waters, and that degassing is higher during high water table periods.Indeed, degassing in streams is higher during high flow periods (0.90 mmol m -2 d -1 ) than base flow periods (0.40 mmol m -2 d -1 ) as a consequence of both higher inputs of groundwater DIC in streams and higher water turbulence.As a matter of fact, degassing is a function of river flow that induces water turbulence and thus increases the gas transfer velocity (Zappa et al., 2007;Raymond et al., 2012).Finally, during the sampling period mean degassing represented 30 % of the total carbon export (i.e., and 60 % of the DIC export) and thus, a significant part of the carbon exported from forest plot rapidly returns in the atmosphere in the form of CO 2 through degassing process, while export of total carbon represents 1.5 % of the mean NEE.* Climate change may affect precipitation amount and its temporal distribution and thus groundwater recharge in ecosystems.
Our study suggests that drought period might enhance soil respiration (heterotrophic + autotrophic) and reduce forest productivity, leading to higher groundwater CO 2 concentration and hence higher degassing of CO 2 in streams.We believe that the probable increase of CO 2 degassing in streams associated with the decrease of forest productivity caused by increasing drought period is a possible transitory positive feedback on atmospheric carbon concentration and has to be taken into account in climate model.As suggested by Loustau et al. (2005), where detrimental effects on forest productivity are expected, enhancement of drought resistance (e.g., through species substitution) may limit the restriction to forest growth.

Conclusion
The work presented here extends the scope on the interaction between the net productivity of forest ecosystem and the concentrations of dissolved carbon in the underlying groundwater.Using net ecosystem exchange measurements, dissolved carbon concentrations determination in both groundwater and stream water and hydrological measurements, we have shown that main controls of carbon fluxes are from hydrology, in interaction with ecophysiological activity of plants (Tab.5, Fig. 7).Thus, high water table resulted in larger DIC and DOC exports and the height of water table outcome from the balance between precipitation, evapotranspiration, water storage and drainage, with a notable influence of local topographic heterogeneities.This work also revealed a large difference of dissolved carbon concentrations between groundwaters and stream waters, those differences originates from distinct processes (e.g., respiration of DOC in groundwater, degassing of CO 2 in streams).Finally, we found that mean export of total dissolved carbon through groundwater represents only 1.5 % of the mean NEE, while 30 % of this carbon export returns to the atmosphere in the form of CO 2 through degassing in first order stream waters.
m are the drainage of groundwater and the mean concentration of DIC in groundwater between two sampling days, respectively in m d -1 and mmol m -3 .Biogeosciences Discuss., doi:10.5194/bg-2017-90,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.We calculated DOC ex as the same manner as DIC ex .Finally we calculated the residence time of DIC (DIC rs ) in groundwater relative to their mean stock (S DIC ) as: DIC rs = S DIC / outputs fluxes (6) Where, Outputs fluxes are DIC ex plus DIC s (only when DIC s is negative) and expressed in mmol m -2 d -1 and, -2015) extremely high precipitations occurred in Jan. 2014 (247 mm compared with an average rainfall of 77 mm for the 2000-2015 period), Sep.2014 (5.5 mm compared with an average rainfall of 69 mm for the 2000-2015 period) and Dec. 2015 (5.7 mm compared with an average rainfall of 95 mm for the 2000-2015 period) received conversely an extremely low amount of precipitations (Fig 2).
g C m -2 d -1 that is equal to 160 mmol m -2 d -1 ) or more recently byHotchkiss et al. (2015) (i.e., 0.87 g C m -2 d -1 that is equal to 72.5 mmol m -2 d - 1 ).Craft et al. (2002) also determined respiration rates value (range 6.4-210 mmol m -3 d -1 ) within a floodplain aquifer of a large gravel-bed river in north-western Montana (USA).Between 12/02/2014 and 16/05/2014 the mean height of the saturated zone is 9.5 m at the Bilos site that leads to a respiration rate of 9.8 mmol m -3 d -1 , consistent with findings ofCraft et al. (2002).Contrary to deeper soil, dissolved organic carbon in the upper soil horizons generally consists in labile low Biogeosciences Discuss., doi:10.5194/bg-2017-90,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.
Biogeosciences Discuss., doi:10.5194/bg-2017-90,2017 Manuscript under review for journal Biogeosciences Discussion started: 12 April 2017 c Author(s) 2017.CC-BY 3.0 License.processes that change DIC or DOC stocks in Bilos groundwater.Overall, residence time (DIC rs or DOC rs ) between two sampling dates that "tends toward zero" indicates negative storage of carbon in groundwater and thus a loss of carbon in groundwater (Fig 5b; 6).In contrast, residence time (DIC rs or DOC rs ) between two sampling dates that "tends toward infinity" indicates positive storage of carbon in groundwater and thus a gain of carbon in groundwater (Fig 5b; 6).

Figure 1 :
Abbreviations Definitions P Precipitation (mm d -1 ) GWS Groundwater storage (mm d -1 ) ETR Evapotranspiration (mm d -1 ) D Drainage (mm d -1 ) H m Mean Bilos Groundwater table (mm) TA Total Alkalinity (µmol L -1 ) pCO 2 Partial pressure of CO 2 (ppmv) DIC Dissolved inorganic carbon DOC Dissolved organic carbon DIC m Mean concentration of DIC in Bilos groundwater (mmol m -3 ) DOC m Mean concentration of DOC in Bilos groundwater (mmo m -3 ) DIC m1 Mean concentration of DIC in first order streams (mmol m -3 ) DOC m1 Mean concentration of DOC in first order streams (mmol m -3 ) DIC s Storage of DIC in groundwater (mmol m -2 d -1 ) DOC s Storage of DOC in groundwater (mmol m -2 d -1 ) DIC ex Export of DIC through drainage of groundwater (mmol m -2 d -1 ) DOC ex Export of DOC through drainage of groundwater (mmol m -2 d -1 ) F Degass F c

Figure 2 :
Figure 2: Seasonal variations of hydrological parameters in the Leyre watershed.Top panels represent the mean monthly variations (temporal variability) of Leyre River flow (left) and Precipitation (right) over the sampling period (2014-2015) and a larger period.(a) Discharge of the Leyre River associated with groundwater table at the Bilos site.(b) Monthly precipitation (P), drainage (D), evapotranspiration (ETR) and groundwater storage (GWS) at the Bilos site.The gray side bars represent high flow periods.

Figure 3 :
Figure 3: Monthly water mass balance (see section 2.5) at the Bilos site for 2014-2015.Pearson coefficient R = 0.85, pvalue < 0.001.Blue points represent months where GWS (Mar.2014, Apr.2014, Mar.2015, Apr.2015, Jun.2015, Jul.2015) is extremely negative (see Fig.2).These blue points are further away from the 1:1 Line than the other months (represented in black).The drainage of the Leyre River is delayed compared to the drainage of Bilos plot.Thus, when the loss of groundwater is extremely high (GWS negative), calculated D do not correspond to measured GWS.Hence, we expected more mistakes when GWS is extremely negative.

Figure 4 :
Figure 4: (a) Mean concentration of DIC (DIC m ) and DOC (DOC m ) in groundwater as a function of water table (H m ).Temporal variations throughout the sampling period of (b) DIC in groundwater (Bilos and the two other piezometers) and DIC in first order streams (medium dashed line) and (c) DOC in groundwater (Bilos and the two other piezometers) and DOC in first order streams (medium dashed line).The gray side bars represent high flow periods.

Figure 5 :
Figure 5: Temporal variations throughout sampling period of (a) ecological parameters at the Bilos site (GPP, R and NEE), here GPP is represented negative (b) storage of DIC and DOC in groundwater at the Bilos site and (c) export of DIC and DOC throughout Bilos groundwater and degassing of CO 2 in first order streams.

Figure 6 :
Figure 6: Residence time of DIC and DOC in Bilos groundwater relative to export and loss (storage when negative).Residence time is calculated as the stocks of each C species divided by their outputs (storage when negative plus export).

Figure 7 :
Figure 7: Relationships between the studied parameters in the Leyre watershed."sign +" represents a positive relationship whereas "sign -" represents a negative relationship.Dashed lines represent the indirect relationship of the water mass balance (P, ETR, GWS and D) on the water table level.

Figure 4
Figure 4 Figure 5 776 Figure 6 778 779 780 781 3 in Sep.2014 and Oct. 2014 respectively) and lowest values of DOC in Bilos groundwater occurred during periods of base flow and low water table (Fig.4).In addition, highest values of DIC (5,400; 5,100 and 3,975 mmol m -3 in Sept, Oct and Nov.2014 respectively) in groundwater were also associated both with highest values of DIC in first order streams, although with lower values (1,300 mmol m -3 in Sep.2014).We also noticed that highest DIC concentrations in the groundwater exactly coincided with positive NEE (R eco > GPP, 21, 33 and 50 mmol m -2 d -1 , respectively in Sep, Oct and Nov 2014) (Fig.4b; 5a).
reported DOC concentrations in soil waters of forest top layers to range from 1,600 to 7,500 mmol m -3 (20 to 90 mg L -1 ).Additionally, in forested very acidic soil worldwide, Camino-Serrano et al. (2014) found DOC concentrations in soil water close to 40 mg L -1 .Hence, combining all these factors, during high groundwater table periods and high flow periods, DOC in Leyre groundwater are very high because of the direct leaching in groundwater of soil organic matter of upper organic horizons (i.e., we found groundwater DOC close to the mean value ofMichalzik et al.