Saturated areas through the lens: 1. Spatio‐temporal variability of surface saturation documented through thermal infrared imagery

Surface saturated areas are key features in generating run‐off. A detailed characterization of the expansion and contraction of surface saturation in riparian zones and its connectivity to the stream is fundamental to improve our understanding of the spatial and temporal variability of streamflow generation processes. In this first contribution of a series of two papers, we used ground‐based thermal infrared imagery for characterizing riparian surface saturation seasonal dynamics of seven different riparian areas in the Weierbach catchment (0.42 km2), a small forested catchment in Luxembourg. We collected biweekly panoramic images of the seven areas over a period of 2 years. We identified the extension of saturation in each collected panoramic image (i.e., percentage of pixels corresponding to saturated surfaces in each riparian area) to generate time series of surface saturation. Riparian surface saturation in all areas was seasonally variable, and its dynamics were in accordance with lower hillslope groundwater level fluctuations. Surface saturation in the different areas related to catchment outlet discharge through power law relationships. Differences in these relationships for different areas could be associated with the location of the areas along the stream network and to a possible influence of local riparian morphology on the development of surface saturation, suggesting a certain degree of intracatchment heterogeneity. This study paves the way for a subsequent investigation of the spatio‐temporal variability of streamflow generation in the catchment, presented in our second contribution.

tial and temporal variability of streamflow generation processes. In this first contribution of a series of two papers, we used ground-based thermal infrared imagery for characterizing riparian surface saturation seasonal dynamics of seven different riparian areas in the Weierbach catchment (0.42 km 2 ), a small forested catchment in Luxembourg. We collected biweekly panoramic images of the seven areas over a period of 2 years. We identified the extension of saturation in each collected panoramic image (i.e., percentage of pixels corresponding to saturated surfaces in each riparian area) to generate time series of surface saturation. Riparian surface saturation in all areas was seasonally variable, and its dynamics were in accordance with lower hillslope groundwater level fluctuations. Surface saturation in the different areas related to catchment outlet discharge through power law relationships. Differences in these relationships for different areas could be associated with the location of the areas along the stream network and to a possible influence of local riparian morphology on the development of surface saturation, suggesting a certain degree of intracatchment heterogeneity. This study paves the way for a subsequent investigation of the spatio-temporal variability of streamflow generation in the catchment, presented in our second contribution.

K E Y W O R D S
catchment hydrology, ground-based thermal infrared imagery, intracatchment variability, remote sensing, riparian processes, surface saturation dynamics, surface saturation mapping 1 | INTRODUCTION Saturation-excess overland flow and its connection to saturated areas were first documented in the seminal work by Dunne and Black (1970). Surface saturated areas (i.e., areas presenting water at the ground surface) have been recognized as key areas for mediating the onset and offset of hydrological connectivity between hillslopes and streams in humid temperate environments (Ambroise, 2004;Bracken & Croke, 2007;Hewlett, 1961;Tetzlaff et al., 2007). Saturation dynamics and the associated hydrological connectivity have been on the agenda of both modelling and experimental studies. Many monitoring studies of surface saturation focused on the near-stream area (i.e., riparian area), which is particularly relevant for run-off generation due to its intrinsic proximity to the stream. Several studies have shown that the spatial extent of near-stream surface saturated areas is a valuable indicator of the general hydrological state of the catchment and, in particular, of groundwater storage during baseflow conditions (i.e., Gburek & Sharpley, 1998;Myrabø, 1997). During precipitation events, riparian surface saturated areas can quickly extend and convey event water to the stream and act as mixing areas for hillslope water contributions (Soulsby, Birkel, & Tetzlaff, 2016;Tetzlaff, Birkel, Dick, Geris, & Soulsby, 2014;Weill et al., 2013). This is often observed in catchments characterized by confined valley bottoms, where persistent saturation can develop in riparian locations with low relief and a shallow water table Niedda & Pirastru, 2014). Thus, an accurate characterization of expansion and contraction dynamics of riparian surface saturation is needed to fully interpret the hydrological behaviour of catchments exhibiting these features and to accurately predict run-off dynamics and associated water and nutrient flowpaths.
Varying riparian morphological traits and upland topographic characteristics have been associated to the variability in hydrological and biogeochemical functions of riparian zones. Several studies have accounted for the spatial variability of riparian zones when exploring riparian functions such as water table fluctuation (Grabs, Bishop, Laudon, Lyon, & Seibert, 2012), vertical and lateral connectivity (Leach et al., 2017;Ploum, Leach, Kuglerová, & Laudon, 2018), or water travel distance and retention of chemicals (Grabs et al., 2012;Ledesma et al., 2018;Vidon & Hill, 2004). On the other hand, riparian surface saturation dynamics have been mainly investigated by taking into account single riparian sections (Zillgens, Merz, Kirnbauer, & Tilch, 2007) or the dynamics of the riparian system as a whole (Ocampo, Oldham, Sivapalan, & Turner, 2006). Due to the spatial variability of riparian characteristics (i.e., riparian width, slope, and soil depth), monitoring of surface saturation that is restricted to a single riparian section can be far from being representative of the whole catchment's riparian zone. Similarly, mapping the dynamics of surface saturation of the riparian zone as a whole may conceal important small-scale variability. Therefore, it is fundamental to characterize riparian surface saturation by accounting for riparian zones' spatial heterogeneity.
As early as in 1975, Dunne et al. made a call for the development of a routine method for the "recognition and quantification of the seasonal and in-storm [and inter-storm] variation of the saturated runoffproducing zones." Progress towards a routine method for mapping the spatio-temporal variability of saturated areas in humid environments remains hampered by technological limitations, especially when it comes to mapping surface saturation dynamics with high spatial and temporal resolutions. In order to get a better understanding of the spatial and temporal scales at which previous studies have addressed surface saturation, we reviewed a total of 64 studies on surface saturation. In 25 of the reviewed studies, surface saturation dynamics were estimated through the use of proxies for surface saturation such as riparian water table level variations (e.g., Waddington et al., 1993;Vidon & Hill, 2004;Ocampo et al., 2006;Tetzlaff et al., 2014), modelling approaches (e.g., Appels, Bogaart, & van der Zee, 2016;Beven & Kirkby, 1979;Blumstock, Tetzlaff, Dick, Nuetzmann, & Soulsby, 2016;Dick, Tetzlaff, Birkel, & Soulsby, 2015;O'Loughlin, 1987;Weill et al., 2013), or a combination of the two (e.g., Baker, Wiley, & Seelbach, 2001;Frei, Lischeid, & Fleckenstein, 2010;Myrabø, 1997;Stieglitz, 2003;von Freyberg, Radny, Gall, & Schirmer, 2014). Proxies for surface saturation such as water table dynamics can be collected at high temporal resolution but are limited to punctual spatial observations.
Modelling approaches for estimating surface saturation dynamics commonly rely on the estimation of topography-based wetness such as the topographic wetness index employed in TOPMODEL (Beven & Kirkby, 1979), multiple existing topographic wetness index variants, or geomorphic indices . These models allow for an estimation of surface saturation over large spatial extends (up to hundreds of km 2 ). However, some of the models' underlying assumptions may not always be valid (e.g., the local slope may not be a valid proxy of the downslope hydraulic gradient), especially in catchments of flat terrain Rodhe & Seibert, 1999). Other modelling studies relied on spatially distributed, physically based simulations of surface saturation dynamics (e.g., Frei et al., 2010;Weill et al., 2013), yet these studies lacked a detailed assessment of the validity of the model results against field observations. Direct mapping of surface saturation (rather than relying on the use of proxies or modelling) was performed within 30 of the 64 reviewed studies for varying spatial extents and with varying monitoring frequencies (Figure 1 provides a graphical representation of the space-time sampling characteristics of surface saturation mapping from these 30 studies, plus this contribution). Except for a few exceptions (i.e., , field surveys (e.g., squishy boot method) have been applied for mapping areas below 5 km 2 and at low temporal resolution (i.e., mainly monthly and punctual observations).
Saturation mapping via remote sensing tools has been mainly relying on satellite and airborne platforms. These techniques are less labourintensive compared with field surveys and can deliver a higher amount of observations in a certain time frame and for larger areas (i.e., > 5 km 2 ; de Alwis, Easton, Dahlke, Philpot, & Steenhuis, 2007;. However, similar to field surveys, remote sensing observations from satellite platforms do not provide the necessary spatial and temporal resolution for detecting heterogeneous riparian surface saturation dynamics within a catchment. Ground-based remote sensing techniques (i.e., thermal infrared [TIR] or visible light imagery) can provide observations at higher temporal (i.e., minutes to weeks) and spatial (i.e., centimetres to meters) resolutions (Glaser et al., 2016;. These techniques will likely become pivotal in generating new, more detailed insights into the functioning of surface saturated area variability and dynamics. Similarly to groundbased TIR, other techniques based on temperature detection (such as thermal imagery from unmanned aerial vehicles and fibre optic distribute temperature sensing) can also provide observation at high spatial (i.e., centimetres to kilometres) and temporal (i.e., minutes to weeks) resolutions, although until today, they have been primarily employed for the characterization of longitudinal stream temperatures and detection of GW exfiltration (Briggs, Dawson, Holmquist-Johnson, Williams, & Lane, 2019;Selker, van de Giesen, Westhoff, Luxemburg, & Parlange, 2006). Within the 64 reviewed studies, 11 did not report clear information on the spatial and temporal scales at which surface saturation was addressed.
Here, we analyse the temporal variability of different riparian surface saturated areas under a new resolution and perspective-namely, "through the lens" of a TIR camera. We employed ground-based TIR imagery in an analogous approach to Glaser et al. (2018; i.e., to detect temperature differences between the water at the ground surface-saturated areas-and the surrounding environment-unsaturated areas), to obtain a unique dataset of 2 years of biweekly observations of different riparian surface saturated areas within the Weierbach catchment in Luxembourg. This long-term studied headwater catchment (0.42 km 2 ) is a reference site for rainfall-dominated mountainous catchments (Zuecco, Penna, & Borga, 2018). The Weierbach is characterized by homogeneous pedology and geology and exhibits a hydrological response that is highly influenced by the wetness state of the system. This leads to a singlepeak response during dry conditions and a double-peak response during wet conditions-after a threshold in catchment storage is exceeded (Martínez-Carreras et al., 2016). The Weierbach catchment has a well-developed riparian zone, characterized by variable morphology (e.g., riparian width and elevation) and the presence of perennial and/or temporary groundwater exfiltration points. Although there is a reasonable understanding of how the overall hydrological response of the catchment is generated Glaser et al., 2016;Klaus, Wetzel, Martínez-Carreras, Ector, & Pfister, 2015 F I G U R E 1 Space-time sampling characteristics of surface saturation mapping from 30 different studies that employed direct mapping of surface saturation, plus this contribution. The size of the studied area refers to the overall area of the investigated catchment or hillslope and the quantity value refers to the total number of times the area has been mapped (information acquired but not directly used in the publication has been included). Studies where different surface saturation mapping methods have been employed or where the same method was employed for different areas have been considered as multiple examples. Methods are indicated with "survey" (e.g., squishy boots method), "remote sensing" (e.g., ground-based and satellite), and "pedo-geo-botan" (i.e., pedological, geological, and botanical aspects used to delineate permanently surface saturated areas of the catchment). The total duration of the mapping period is indicated close to the circles for nonpoint observations (Y = years, M = months, W = weeks). In order to make the studies comparable in terms of surface saturation frequency of observation, we considered only the most recurrent time interval between two observations for each study. Studies that reported a seasonal mapping were included under the monthly frequency. Note that the big circle (which corresponds to a year of digital images acquired every minute during day time by  is not in scale for display purposes. References for the 30 studies considered for the review figure: Bari, Smettem, & Sivapalan, 2005;Blazkova, Beven, & Kulasova, 2002;Brun et al., 1990;Buttle & Sami, 1992;Chabot & Bird, 2014;Creed, Sanford, Beall, Molot, & Dillon, 2003;D. A. De Alwis, Easton, Dahlke, Philpot, & Steenhuis, 2007;Devito, Creed, & Fraser, 2005;Franks, Gineste, Beven, & Merot, 1998;Gineste, Puech, & Mérot, 1998;Glaser et al., 2016;A. Güntner, Uhlenbrook, Seibert, & Leibundgut, 1999;Andreas Güntner et al., 2004;Inamdar & Mitchell, 2007;Kulasova, Beven, Blazkova, Rezacova, & Cajthaml, 2014;McDonnell & Taylor, 1987;Rinderer, Kollegger, Fischer, Stähli, & Seibert, 2012;. The bibliography for the 30 studies considered for the review figure can be found in Appendix A1 et al., 2015), there is still a lack of understanding of the dynamics of small-scale riparian processes, like the spatial and temporal variability of riparian surface saturation, and of how these dynamics are related to the hydrological response (Scaini et al., 2017).
In this first contribution of a series of two papers, we apply ground-based TIR imagery as a routine method for mapping surface saturation dynamics across multiple seasonal and hydrological conditions and across multiple sites in the riparian zone of the Weierbach catchment. We analyse the spatio-temporal dynamics of surface saturation by applying statistical analyses on an extensive surface saturation dataset produced by direct field observation. In particular, through this novel approach, we investigate the following questions on saturated area dynamics: 1. Are the overall surface saturation dynamics (i.e., seasonal and yearly dynamics) of the seven investigated areas similar? 2. How do hydrological conditions (i.e., precipitation, stream discharge, evapotranspiration, groundwater level, soil moisture, and catchment storage) relate to the temporal (seasonal and yearly) variability of surface saturation in different riparian locations in a catchment?
We leverage the outcomes of this study in the accompanying manuscript for investigating how hillslope-riparian-stream (HRS) connectivity is established in the Weierbach catchment. This will eventually improve our understanding of how the spatial variability of streamflow generation is linked to surface saturation dynamics.
Elevation ranges from 458 to 513 m.a.s.l. Topography is characterized by a quasihorizontal plateau, covering 54% of the catchment and cut by steep (≥5 ) V-shaped valleys.
Vegetation is composed of Oak and Beech trees on the western side of the catchment and Spruce on the eastern side. Ferns and herbaceous plants dominate in the riparian zone. In this study, we identify the riparian zone considering a combination of different criteria.
The change in dominant vegetation and the presence of shallow clayloam, organic soil (i.e., Leptosol), peculiar of the low relief near-stream area of the catchment, set a visual basis for differentiating riparian from other landscape elements (i.e., hillslopes, plateau). The riparian zone is gently sloped (<5 ) and covers 1.2% of total catchment area.
The catchment's run-off response to precipitation is influenced by a storage threshold (Martínez-Carreras et al., 2016) and changes between the dry and the wet seasons. In case of dry antecedent conditions, the catchment produces a single spiky peak of short duration (i.e., hours), whereas the response is bimodal during wet antecedent conditions-with a first peak followed by a broader second peak of longer duration (extending up to several days). Martínez-Carreras et al. (2016) showed that the first peak is mainly composed of water from precipitation, throughfall, and rapid HRS connectivity through F I G U R E 2 Location and instrumentation map of the Weierbach catchment saturation excess overland flow and preferential flowpaths (such as macropores and/or fractures along the hillslopes) whereas the second peak mainly consists of infiltrated soil water and groundwater flowing though the fractured bedrock, once the storage threshold is exceeded.

| Hydrometeorological measurements and catchment storage calculation
Hydrometeorological measurements are carried out in the Weierbach catchment since 2002 as part of a long-term monitoring programme.
Water levels were measured and recorded every 15 min by a pressure transducer (ISCO 4120 Flow Logger) installed at a V-notch weir at the catchment outlet ( Figure 2) and translated into discharge via a rating curve (based on salt dilution measurements). Precipitation data were measured at a canopy-free location in the Weierbach catchment (see Catchment storage estimates were calculated following the methodology developed by Martínez-Carreras et al. (2016). In their study, the total amount of water stored in the catchment at a given time was calculated as the sum of storage in three separate zones: where S UNSAT is the water stored in the variably unsaturated zone (estimated from the VWC of the soil above the water content at field capacity-measured at the different soil depths in Sites 3, 4, 5, and 7), S SAT is the water stored in the variably saturated zone (within the range of water table fluctuations-estimated from GW levels in the hillslope and in the plateau: GW3 and GW5, respectively, in this study), and S RES is the water stored in the residual saturated zone (i.e., estimated drainage porosity of the basal layer, fractured bedrock, and fresh bedrock). S TOTAL is obtained for the different landscape elements by multiplying the value of total storage by the area of each element. For more specific information on the calculation of catchment storage, the reader is referred to Martínez-Carreras et al. (2016).

| Monitoring of surface saturated areas in the riparian zone
We focused on seven distinct riparian areas in the Weierbach catchment ( Figure 3). Each area was labelled with an abbreviation (cf. areas' name in Figure 3) indicating the stream branch where it is located (i.e., L = riparian areas on the left stream branch; M = riparian areas on the middle branch; R = riparian areas on the right branch, and S = riparian areas on the main stream) and its position along the branch (i.e., numbered from downstream to upstream). Descriptive topographic characteristics of the different riparian areas, such as average elevation, area extension (i.e., area covered by the monitoring), and maximum riparian width, were extracted from a high resolution LIDAR DEM (~5-cm resolution). TIR observations (i.e., sequential images and videos) and visible light photographs of the different riparian areas were acquired for a total of 63 mapping campaigns with a weekly to fortnightly recurrence interval from November 2015 to December 2017. The used handheld TIR camera (FLIR T640, FLIR Systems, Wilsonville, OR, USA) is sensitive to the radiation emitted from an observed surface (or the first 0.1 mm of a water column) over a spectral range of 7.5 to 14 μm, produces images of 640 × 480 pixels, and covers a temperature range of −40 C to 2000 C, with a thermal sensitivity of <0.035 C at 30 C. Information about object emissivity "ε" (usually set between 0.95 and 0.97 for freshwater), atmospheric temperature, air humidity, object's distance from the device, and reflected ambient temperature were provided to the camera in order to correct the detected temperature for these parameters (the correction is automatically done by the camera's software).
The final product of the TIR camera is an image (or video) reporting surface temperatures for each image pixel. This temperature information can be used to classify the pixels into pixels corresponding to water ponding or flowing at the ground surface (saturated pixels, i.e., stream and riparian ponds) and pixels representing surrounding material (unsaturated pixels, i.e., soil, rock, and vegetation; Figure 4a,b). In order to be able to discern these two pixel classes, a clear temperature contrast between surface water (saturated pixels) and surrounding material (unsaturated pixels) is required.
Moreover, the camera view on the saturated surfaces has to be free from obstructions (e.g., vegetation, snow, fog, and heavy rain). Below, we shortly explain how we transformed the information from the acquired TIR images into information on the extent of surface saturation. A detailed description of the postprocessing workflow and the TIR imagery technique applied for surface saturation mapping in general can be found in Glaser et al. (2018).
In order to prepare the TIR images (or videos) for the extraction of the extent of surface saturation in the investigated riparian areas, we followed a sequence of postprocessing steps. The sequence consisted of (a) creating panoramic images by overlapping single images (or video frames; Figure 4a Glaser et al. (2018), where the methodology for the TIR approach was developed. We then calculated the percentage of saturated pixels in each panorama as proxy for the extent of surface saturation in the investigated areas following the manual approach for the generation of saturation maps described in Glaser et al. (2018). This approach consists in (a) manually selecting the temperature range corresponding to surface saturation, (b) adapting the selected range to create a saturation map with a pattern of saturated pixels matching best the saturation pattern identified via visual inspection of the TIR panoramas and visible light images (here defined as optimal solution), and (c) calculating the number of pixels falling into that temperature range over the total number of pixels of each image (Figure 4e). The three steps were repeated for each location and observation date, meaning that an individual temperature range was selected for each TIR panorama. TIR panoramas that showed poor temperature contrast and/or high influence from obstructing elements were excluded from the analyses (34% of the 441 acquired panoramas).
As shown by Glaser et al. (2018), the manual selection of temperature ranges is to date the best approach for generating reliable saturation maps from TIR datasets where images show very variable conditions (e.g., in terms of wetness and overall temperature range) and present slight perspective shifts. However, because the manual selection of an optimal solution for the saturation estimation is a subjective process, different operators may tend to select different optimal temperature ranges, including more or less pixels into the group of saturated pixels based on their individual perception.
F I G U R E 3 Location and example of visible light and thermal infrared (TIR) panoramic photo of each investigated riparian area for a wet and dry condition. Visible light images are shown for wet conditions only. In the TIR panoramas collected during wet conditions, saturated pixels correspond to the lighter colours. In the TIR panoramas collected during dry conditions, saturated pixels correspond to the darker colours (except for Area L1; pictures: M. Antonelli and B. Glaser). Yellow arrow: flow direction. Red circles: location of permanent springs To investigate the range of possible surface saturation outcomes, we varied for some panoramas the width of the saturated pixels temperature range (i.e., changing the higher and the lower values in small temperature steps) until the saturation pattern clearly mismatched the saturation pattern selected as the optimal solution (cf. Glaser et al., 2018). The saturation pattern including the higher number of pixels, and still reflecting the realistic saturation pattern, was used to determine the maximum estimate of surface saturation. The saturation pattern including the lower number of pixels, and still reflecting the realistic saturation pattern, was used to determine the minimum estimate of surface saturation. We estimated different saturation outcomes for the investigated riparian areas taking into account images collected during different saturation levels (i.e., five to seven images for each investigated area). We then plotted the minimum and maximum estimates of saturated pixels against the optimal estimation (within each area) and determined a regression equation from which we retrieved the minimum and maximum estimates of saturation for the whole time series of saturation of each of the areas.
The overall amount of surface saturation estimated from the TIR images in each area represents both riparian surface saturation and water in the stream channel. The stream channel receives water contributions from the riparian zone (i.e., lateral contribution; cf. Figure 3 red circles) and water supply from upstream along the stream channel itself (i.e., longitudinal/upstream contribution). The relative amount of water provided by these two different contributions is difficult to disentangle. However, between the investigated riparian areas, we expect upstream contributions to be higher in downstream areas compared with headwater areas. In order to strengthen the comparison of the relationships between riparian surface saturation in different areas (i.e., headwater areas vs. downstream areas) and stream discharge (cf. Figure 11 in Section 4), we provide an example of an estimation of downstream areas' surface saturation with a reduced influence of the F I G U R E 4 Workflow of the thermal infrared (TIR) image postprocessing for the example of the TIR panorama of Area M2 of February 25, 2016. For more details, the reader is referred to Glaser et al. (2018) upstream contribution (i.e., excluding the stream from the calculation of surface saturation). However, this exercise leads to a considerable loss of information on the overall level of surface saturation because, by doing so, a substantial part of surface water represented by the lateral inflows is also excluded. For this reason, the estimation of surface saturation with a reduced influence of the upstream contribution is only provided as an example, and it is not employed in all the analyses.

| Statistical data analysis
We analysed the time series of saturation of the seven riparian areas, in order to investigate the temporal dynamics of surface saturation.
We applied a min-max normalization to the percentage of saturated pixels for each area. We expressed the normalized values as a percentage in order to compare areas of different extension (cf. Section 4.1). We will refer to these values as "normalized saturation." For the normalization, we accounted for the percentage of saturated pixels from images acquired during periods where surface saturation was not affected by particular meteorological conditions such as frozen soils (which will be represented by normalized percentages below 0%) or rain-on-snow events (which will be represented by normalized percentages above 100%). By quantifying the observations obtained during the occurrence of frozen soils and rain-on-snow events in this way, they can be easily identified in the figures and provide information on the field conditions. These observations where retained in the statistical analysis of the dataset because they are not statistical outliers. Nevertheless, we tested the statistics excluding F I G U R E 5 Time series of (a) precipitation (black) and reference evapotranspiration (grey-smoothed trend of reference ET showed in red), (b) discharge and catchment storage (red arrows show moments when discharge response was more pronounced than catchment storage), (c) soil volumetric water content along the hillslope-riparianstream transect (10-cm depth), and (d) ground water levels from October 2015 to January 2018 these observations and the results remained consistent. Descriptive statistics and the nonparametric Mann-Whitney-Wilcoxon Test (α = .05) were used to compare the temporal distribution of the normalized saturation in different areas.
In order to explore a possible influence of precipitation and evapotranspiration on the seasonal dynamics of surface saturation in the different areas, we compared the double mass curves (DMCs) of rainfall and run-off of the catchment for the two investigated hydrological years (HYs) with the DMCs of rainfall and surface saturation in the different areas (i.e., we cumulated the estimated values of normalized surface saturation from one date of observation to the other and highlighted moments of vegetation growth and high evapotranspiration).
Classical rainfall-run-off DMCs can provide direct information on seasonal run-off formation (Pfister, Iffly, Hoffmann, & Humbert, 2002;Seibert, Jackisch, Ehret, Pfister, & Zehe, 2017). Similarly, by evaluating how cumulated surface saturation evolves in response to cumulated rainfall, we aimed to obtain information on the influence of seasonal variables (i.e., precipitation and evapotranspiration) on the development of surface saturation. Note that by cumulating normalized surface saturation, we do not intend to give an estimation of a total amount of surface saturation of each HY. Instead, we consider the cumulated surface saturation as a way to identify periods of general increase or decrease of saturation.
As a measure of how fast the saturation changed in each area, we calculated the difference between the normalized saturation estimated on two consecutive dates and divided this value by the number of days in each period to obtain daily normalized rates of change. We tested similarities between the daily rates of change between different areas with the Mann-Whitney-Wilcoxon Test (α = .05). The test was applied by taking into account each surface saturated area against every other area for the dates when an estimation of a change rate of saturation was available for both areas. Additionally, we tested if the difference in normalized saturation estimated between two consecutive dates was related to differences in GW levels, soil VWC along the HRS transect and soil VWC profiles at Sites 3, 4, 5, and 7, catchment storage or to the amount of precipitation (expressed via the antecedent precipitation index-as per McDonnell, Owens, & Stewart, 1991) observed between the same dates. We applied Spearman's rank F I G U R E 6 Time series of soil volumetric water content measured at 10-, 20-, 40-, and 60-cm depth at Sites 3, 4, 5, and 7, from October 2015 to January 2018 correlation to test these relationships (α = .01). As before, we only took into account the dates for which an estimation of saturation was available for the analysed area.
We applied Spearman's rank correlation test rho (ρ; α = .01) in order to test monotonic relationships between (a) the time series of normalized saturation estimated in the different investigated riparian areas and (b) between these values and the time series of hydrometric measurements (i.e., daily-averaged values of outlet discharge, estimated catchment storage, GW levels, soil VWC along the HRS transect, and soil VWC profiles at Sites 3, 4, 5, and 7). These relationships were tested for the whole study period.
In order to analyse the shape of the relationship between the surface saturation in the different areas and baseflow discharge at the catchment outlet, we relied on the observations of surface saturation that were not impacted by the occurrence of precipitation (i.e., images taken while rainfall occurred, during rising limbs or peaks of discharge, and at the early stage of discharge recession) or by the occurrence of particular meteorological conditions such as frozen soils or rain-onsnow events. This set of data describes the evolution of surface saturation along the gradual change in wetness state of the catchment and can be related to the surface saturation versus baseflow discharge relationship described by . We fitted various types of equations on the observations not impacted by the occurrence of precipitation, and we found that power law equations (Sat = a*Q b ) adequately approximate the observed trends (fitting carried out on nontransformed data; goodness-of-fit was tested with Kolmogorov-Smirnov test-p value >.1).  Figure A1 T A B L E 1 Characteristics of the investigated riparian areas April 2017) during the HY 2016/2017. In December 2017, a rain-onsnow event produced a high peak discharge. Note that in January 2017, the stream was partially frozen (as a result, no discharge data are available for that period).

| Hydrological response and catchment storage
Shallow soil VWC (10-cm depth) gradually decreased from the riparian zone towards the hillslopes. Riparian soil VWC was oscillating between a maximum of~70% during wet conditions and a minimum of~65% during dry conditions. Shallow soil VWC at the other monitored locations was more variable (Figure 5), showing marked reaction to precipitation. Soil VWC of the shallow soil measured in the Spruce-covered hillslope revealed a tendency of the soil to dry more rapidly than in other locations. Soil VWC measured in Sites 3 and 4 was generally more responsive to precipitation compared with soil VWC measured at Sites 5 and 7, at all depths ( Figure 6). Soil VWC measured at Site 3 decreased form the 10-cm depth to the 60-cm depth, whereas the opposite was observed at Site 7. Soil VWC measured at Sites 4 and 5 was more similar along the depth profile ( Figure 6). As previously observed, soil VWC profiles measures in the Spruce-covered and Pine-covered hillslopes at all depths revealed a tendency of the soil to dry more rapidly than in other locations. GW . Out of the total number of 441 acquired TIR panoramas, 291 panoramas were used for the estimation of a value of surface saturation. Data represented with an asterisk refer to estimated values of surface saturation from TIR panoramas with optimal temperature contrast and no obstructive elements between the camera and the object (n = 101). Data represented with a circle and a triangle refer to TIR observations that were slightly influenced by the presence of vegetation or snow (n = 110) and where the temperature contrast was not optimal (n = 80), respectively, but that were still usable for the estimation of a value of surface saturation. Normalization was done according to the highest and lowest observed percentage of saturation within each area individually, conditions with frost and rain-on-snow excluded. Frost = condition with frozen stream and riparian soil. ROS = rain-on-snow event levels responded to precipitation in all four wells. Water levels in GW2 (close to a source area) and GW3 (hillslope foot) were always shallower than 2.00 m. GW5 and GW6 (both located on plateaus) behaved differently-mostly during dry periods, when the water level in GW5 would recede at a constant rate, whereas GW6 would dry quickly and show a reaction to new water inputs more similar to the response in soil moisture (also in comparison with the shallower GW2 and GW3).
During the wet periods, catchment storage and discharge showed very similar trends. During the dry periods-when catchment storage mainly corresponded to the GW reservoir-precipitation triggered more pronounced changes in discharge (as a single peak) in comparison with storage (cf. Figure 5b, red arrows).

| Characteristics of the investigated riparian areas and correspondent upslope catchments
We assigned the investigated riparian areas in the Weierbach catchment to three main groups, based on intrinsic area characteristics (Table 1). Areas L1, M3, and R3 correspond to the most upstream locations (i.e., source areas) of the stream. They are wide (8.3 m in average) and fed by perennial groundwater exfiltration (e.g., stable exfiltration points observed via TIR imagery throughout the year, cf. Figure 3). Areas M2 and S2 display similar characteristics but are located further downstream (Figure 3). In Areas M1 and R2, the riparian zones are narrower (4.9 m in average) and without clearly identified points of perennial groundwater exfiltration. On the basis of these differences, we qualify the first group as "Stream Source Areas with Perennial springs (PSA),", the second group as "Areas along the stream with Perennial Springs (PSpA)," and the third group as "Areas along the stream with Non-Perennial Springs (N-PSpA)."

| Spatio-temporal dynamics of surface saturation and their relationship with meteorological conditions
The values of the estimated normalized surface saturation were highly monotonically related between the different areas over the whole F I G U R E 9 Violin plots of the distribution of normalized surface saturation from the time series of the seven studied areas (grouped according to the PSA, PSpA, and N-PSpA classification: PSA = Stream Source Areas with Perennial springs; PSpA = Areas along the stream with Perennial Springs; N-PSpA = Areas along the stream with Non-Perennial Springs). Shaded areas highlight values of saturation acquired during particular boundary conditions such as frozen riparian soils (normalized saturation below 0%) and rain-on-snow events (normalized saturation above 100%).  Figure 9). Areas M1, R2, and S2 had particularly high median values (Figure 9; Table 2). The variability of the observations around the mean values (i.e., standard deviation) was similar for all areas (Table 2), although slightly higher in Area S2 (~38%). A summary of the descriptive statistics for the normalized saturation distribution of the seven riparian areas is reported in Table 2

| Relationship between surface saturation dynamics and hydrometric measurements
For all riparian areas, we identified a strong monotonic relationship between normalized saturation and catchment discharge (Spearman's rank test ρ not lower than 0.78 for all areas, p value <.01; Table 3).
We found overall positive and significant monotonic relationships between the normalized saturation in the different areas and GW levels, VWC, and the estimated storage of the catchment (Table 3) Also, VWC measured at Site 4 (low hillslope Spruce covered) was generally less correlated with normalized surface saturation in all investigates riparian areas.
Changes in the extent of surface saturation between two observation dates were significantly related to the changes of GW level in Location GW3 (0.60 < ρ < 0.77). All the areas except Areas M2 and R2 showed also significant correlation with changes of GW level in GW2 (0.58 < ρ < 0.67; Table 4). Area M1 was particularly correlated also with the GW levels measured in Location GW5 (ρ = 0.72).
Changes in the extent of surface saturation between two observation dates were also significantly related to the changes in catchment storage (0.65 < ρ < 0.84) in all areas. A low but significant correlation of the changes of surface saturation with the antecedent precipitation index was observed for the Areas L1, M1, R2, and S2 (ρ = 0.50).
Changes in soil VWC along the HRS transect (10-cm depth) were not significantly correlated to changes in surface saturation between two observation dates for any of the investigated riparian areas (Table 4).
Changes in soil VWC measured at Sites 3, 4, 5, and 7 at 10-, 20-, 40-, and 60-cm depth were significantly correlated to changes in surface saturation (except for Area L1 and soil VWC at Site 3 at 10-cm depth; Table 4). However, only surface saturation changes in Areas R3 and S2 showed a good correlation with changes in soil VWC measured at different sites and depths. In particular, surface saturation changes in Areas R3 and S2 were particularly related to changes in soil VWC measured at Sites 4 and 5 (at low and middle slope positions, respectively) with ρ ≥ 0.56 for Area R3 and ρ ≥ 0.60 for Area S2.
The surface saturation versus outlet baseflow discharge relationship has been investigated for the seven riparian areas (Figure 11). In Areas S2, M1, and R2, the observations obtained during low flow appeared more scattered than for the other areas, probably due to stream water contribution from upstream. In order to give an example of a reduced influence of upstream water contribution in Areas M1, R2, and S2, we re-estimated the extent of surface saturation after having excluded stream pixels from the TIR images. One can notice a general shift towards lower saturation for Areas M1 and R2 (see Figure 11, M1, R2, and S2 "no upstream contr."). This is less pronounced for Area S2, probably as a result of the presence of permanent springs within this area, which maintain the riparian zone generally wetter. In general, scattering in the observations at low flow appears reduced in Areas S2, M1, and R2 after having reduced the influence of upstream water contribution. Note that we did not apply the exercise of reducing upstream contribution to Area M2, because we frequently inferred from the TIR images that the portion of stream T A B L E 3 Spearman's rank correlation between the normalized saturation in the seven areas and hydrometric variables with Area L1 and the other areas in general. When considering the observations affected by precipitation (Figure 11, grey dots), a slight hysteretic effect could be observed for the surface saturation versus discharge relationship, in particular for Areas L1, M2, M3, and R3 at higher discharge.

| DISCUSSION
We have used ground-based TIR imagery for mapping the spatiotemporal dynamics of riparian surface saturation expansion and contraction in the Weierbach catchment. For the first time, the dynamics of surface saturation in different riparian locations within the same catchment have been monitored at a temporal resolution high enough to characterize their seasonal variability. To the best of our knowledge, prior to this study, extensive time series of surface saturation dynamics have been displayed only as model outputs, often considering the overall amount of saturation in the catchment and rarely being validated Weill et al., 2013). We observed strong similarity in the expansion/contraction dynamics between the seven riparian surface saturated areas over the whole study period, although there were some differences in the timing of maximum or minimum levels of saturation in the seven areas. N-PSpA and Area S2 showed generally higher normalized surface saturation values (i.e., high median value). The maximum of surface saturation in these areas occurred during a rain-on-snow event in December 2017. This is likely due to the fact that these areas receive the highest contributions of stream water from upstream than other areas (i.e., Area M2 and PSA; cf. areas' locations in Figure 3). We observed the lowest surface satu- and less water in the stream channel) and the fact that surface water in the riparian zone was frozen. We did not consider frozen surface water in the riparian zone as being surface saturation, because frozen (solid) water exhibits different characteristics than free (liquid) water, for example, in terms of reaction to incident precipitation or movement dynamics.
The yearly and seasonal variability in the dynamics of surface saturation in the seven areas was found to reflect the yearly and seasonal variability of catchment run-off (cf. Figures 8 and 10). Increasing cumulated amounts of precipitation caused cumulated run-off and surface saturation to increase in a similar way during the wet periods.
Increasing ET losses during the vegetative period led to moments of low run-off and low surface saturation (i.e., flatter cumulated run-off and surface saturation) despite that the amount of precipitation did not change considerably. Occurrence of breaks along the DMCs when passing from wet to dry conditions-and vice versa-were very similar between the rainfall-run-off and the rainfall-surface saturation DMCs. However, breaks and slope changes in the rainfall-surface saturation DMCs were generally less sharp than the slope changes observed in the rainfall-run-off DMCs. Martínez-Carreras et al. (2016) showed that the Weierbach catchment's run-off response is influenced by a storage threshold that, once exceeded, allows high T A B L E 4 Spearman's rank correlation between the changes in the amount of surface saturation between two consecutive observation dates and changes of hydrometric variables (GW = ground water; Stor_tot = total catchment storage; Site 3 10 cm = soil volumetric water content measured at Site 3 at 10-cm depth-similar naming for the other sites and depths) and antecedent precipitation index (API) calculated between the same observation dates Note: All shown correlations are significant with α = .01. Changes in soil volumetric water content measured at 10-cm depth along the hillslope-riparian-stream transect did not show significant correlation with changes in surface saturation. discharge volumes to be generated by the catchment even in response to relatively small precipitation events. The similarity between the break points in the run-off and surface saturation DMCs may indicate that the same storage threshold influences the seasonal transition between low and high extents of surface saturation in the riparian areas. However, other aspects may play a role in regulating the seasonal expansion and contraction of surface saturation as well.
For example, the smoother slope changes in the surface saturation DMCs than run-off DMCs may reflect that riparian soil hydraulic characteristics influence the expansion and contraction of surface saturation by defining the degree of resilience of surface saturation to develop in response to increasing and decreasing catchment's wetness conditions. In this sense, the seasonal transition of riparian surface saturation may be subjected to a second, different threshold, which is defined by the riparian soil capacity to store and release water (Zehe, Lee, & Sivapalan, 2006).  (2003) found that the expansion of the saturated area was consistent with GW level dynamics in the lower hillslope and hollow zones in the Maimai catchment. Van Meerveld, Seibert, and Peters (2015) noted in the Panola catchment that the hillslope and the riparian zone only became connected when GW levels rose in the lower part of the hillslope.
We did not observe a consistent relationship between changes in surface saturation and changes in soil VWC measured at different soil profiles. However, changes in surface saturation in Areas R3 and S2 showed good correlation with changes in soil VWC measured in Sites 4 (low hillslope) and 5 (midhillslope) at all depths (Table 4) The saturation-baseflow discharge relationships observed in the different riparian areas can be related to the dynamics illustrated in the perceptual model in Figure 12. At low flow, the differences observed in the saturation-baseflow discharge relationships (i.e., amount of surface saturation and scattering in the observations, cf. Figure 11) can be explained by the presence of perennial springs and the location of the riparian area (i.e., area elevation-which determines the variable amount of water reaching the area from upstream locations). At higher flow, the possibility for saturation to develop upstream in PSA (cf. Figure 12b) could explain the fast change in saturation with increasing baseflow discharge in these areas compared with the others (cf. Figure 11). Indeed, it also has been observed by others that saturation that develops in previously dry channels is more reactive than saturation in riparian areas, which is rather influenced by the speed at which the soil drains . The development of more persistent saturation in the riparian soils than the stream channel (Figure 12c) may explain the slight hysteretic effect that was observed in the saturation-discharge relationship of PSA and Area M2 (cf. Figure 11-grey dots). The hysteretic relationships between saturation and discharge that we observed in some areas provide a first actual feedback to the possible hysteretic relationship between surface saturation and outlet discharge that has been usually observed through modelling approach (Glaser et al., 2016;Weill et al., 2013). However, the hysteretic relationships between saturation and discharge observed in this study were never as clearly defined as those observed in modelling studies (Frei et al., 2010;Weill et al., 2013), despite the relatively high number of observations at high flow stages. In this sense, TIR observations at a higher temporal resolution during precipitation events would help to clarify the hysteretic patterns that may occur in the different areas and could be used to validate hysteretic behaviour observed through simulations. Overall, the small but noticeable differences observed in the saturation-baseflow discharge relationships provided information on the different potential for lateral and longitudinal hydrological connectivity to be established through the different riparian areas during different flow stages.
To date, surface saturation-baseflow discharge relationships have been inferred considering only the total surface saturation extent in a catchment . The surface saturation-baseflow discharge relationship has been defined by Ambroise (2016) as a characteristic curve of the catchment, fundamental for understanding and modelling the interaction of water from different sources on the saturated areas and its influence on streamflow during baseflow conditions. By repeatedly monitoring the dynamics of surface saturation in different areas, we found indication of possible intracatchment variability of this relationship. Moreover, the frequency at which we observed surface saturation in this study F I G U R E 1 2 Proposed perceptual model for the development of surface saturation in PSA (Stream Source Areas with Perennial springs) and N-PSpA (Areas along the stream with Non-Perennial Springs) riparian areas allowed us to explore the dynamics of surface saturation during baseflow conditions under different flow stages. It also allowed us to consider how seasonality may affect the observed dynamics. Considering the fact that the broadly used topography-driven indices and geomorphic indices for estimating surface saturation are known to perform relatively poorly during low flow stages Western, Grayson, Blöschl, Willgoose, & McMahon, 1999), our observation of surface saturation dynamics during low baseflow conditions is particularly valuable for obtaining new insights into riparian processes and potentially improve these indices. In example, from the analysis of the rainfallsurface saturation DMCs, we observed that increasing ET losses during the vegetative period lead to moments of low surface saturation despite the amount of precipitation did not change considerably.
The neglection of this shift in dominant processes in the indices calculation might be the reason for the poor performance during dry periods.

| CONCLUSIONS
This study is a contribution to the call for the development of a routine method for mapping surface saturated areas  and to the need to start characterizing the spatial and temporal variability of riparian processes for a better understanding of catchments hydrological and biochemical functioning (Grabs et al., 2012;Ledesma et al., 2018;Tetzlaff et al., 2008;Vidon & Hill, 2004). We applied TIR technology as a valid routine method for repeated mapping of surface saturation (in our case, at weekly or biweekly frequency) in the Weierbach catchment. The frequency at which we monitored surface saturation was critical to characterize the similarities and differences in both the temporal dynamics of surface saturation in different areas and their relationship with stream baseflow discharge.
The observed yearly and seasonal dynamics of surface saturation in the different riparian areas of the catchment were found to be similar. Based on the analysis of DMCs for the surface saturation in comparison with the DMC of discharge, we hypothesized that storage thresholds control the transition between low extents of surface saturation and high extents of surface saturation in the Weierbach catchment. Another similarity between the dynamics of surface saturation observed in different investigated areas has been found in their relationship with the variability in catchment's storage and GW levels measured in lower hillslope locations. This supports the role of riparian surface saturation as a valuable indicator of groundwater storage during baseflow conditions previously assessed in different studies (i.e., Gburek & Sharpley, 1998;Myrabø, 1997).
The shape of the relationship between surface saturation and baseflow discharge could be approximated with a power law in all cases. However, small differences in the relationships for the different areas could be associated with the location of the areas along the stream network (i.e., area elevation) and with the local riparian morphology (i.e., area width and the presence of GW exfiltration points).
These characteristics represent a source of intracatchment variability that may have implications on the potential of different riparian surface saturated areas in mediating hydrological connectivity along the HRS continuum.
Based on our findings and conclusions, we may now ask "Are all riparian zones in our catchment the same, or would the small differences in their dynamics of surface saturation mirror the degree of hydrological connectivity of the different areas with the hillslopes?" With this question in mind, we will present our investigation on the spatial heterogeneity of streamflow generation in our second contribution.
The data and information obtained in this study will prove essential for investigating the spatial variability of streamflow generation in the Weierbach catchment and its relationship with surface saturation. The same approach used in this study can be potentially employed in other catchments as well, especially in those where the riparian zone represents an important interface between the hillslopes and the stream.

ACKNOWLEDGMENTS
We would like to thank Jean-Francois Iffly, the Observatory for Cli-

CONFLICT OF INTEREST
The authors declare that they have no conflict of interest.

DATA AVAILABILITY STATEMENT
Data used in this study are property of the Luxembourg Institute of Science and Technology. They are available upon request from the authors.