Effect of DEM-smoothing and -aggregation on topographically-based flow directions and catchment boundaries

: It is generally assumed that, in humid climates, the groundwater table is a subdued copy of the surface topography. However, the general groundwater table is unlikely to be affected by the microtopography as seen in high-resolution Digital Elevation Models (DEMs). So far, there has been little guidance on the best resolution DEM to use to determine the shape of the water table or the direction of shallow groundwater flow in headwater catchments. We, therefore, looked at the effects of DEM-smoothing and -aggregation on the calculated flow directions and derived catchment boundaries, and identified areas and landscape features for which the calculated flow directions are particularly sensitive to DEM smoothing or aggregation. For > 40 % of the area of the Krycklan study catchment, the calculated flow directions depend strongly on the degree of smoothing or aggregation of the DEM. The four main landscape features for which DEM smoothing or aggregation strongly affected the calculated flow directions were: local slopes in the opposite direction of the general slope, flat areas, ridges, and incised streambanks. To determine the effects of the changing flow directions on the derived catchment boundaries for the smoothed and aggregated DEMs, we calculated the drainage area for 40 locations, representing the outlets of catchments of varying sizes. The shape and size of the catchments of first-order streams were most affected by the processing of the DEM. These streams were often almost completely smoothed out during the DEM preprocessing steps. These shifts in catchment boundaries and drainage area would have a large effect on the water balance. This study thus highlights the need to carefully consider the effects of DEM smoothing or -aggregation on the calculated flow directions and drainage areas. high-resolution DEM use the shape of the table the of shallow groundwater flow in headwater catchments. at the effects of DEM-smoothing and -aggregation on the calculated flow directions and derived catchment boundaries, and identified areas and landscape features for which the calculated flow directions are particularly sensitive to DEM smoothing or aggregation. For > 40 % of the area of the Krycklan catchment, the calculated flow directions depend strongly on the degree of smoothing or aggregation of the DEM. The four main landscape features for which DEM smoothing or aggregation strongly affected the calculated flow directions were: local slopes in the opposite direction of the general slope, flat areas, ridges, and incised streambanks. To determine the effects of the changing flow directions on the derived catchment boundaries for the smoothed and aggregated DEMs, we calculated the drainage area for 40 locations, representing the outlets of catchments of varying sizes. The shape and size of the catchments of first-order streams were most affected by the processing of the DEM. These streams were often almost completely smoothed out during the DEM preprocessing steps. These shifts in catchment boundaries and drainage area would have a large effect on the water balance. This study thus highlights the need to carefully consider the effects of DEM smoothing or -aggregation on the calculated flow directions and drainage areas.


Introduction
In humid climates, the groundwater table is often assumed to be a subdued replicate of the surface topography (Condon and Maxwell, 2015;Haitjema and Mitchell-Bruker, 2005;Tóth, 1962;Winter, 1999). Digital Elevation Models (DEMs) can, therefore, be used to estimate the direction of groundwater flow (Haitjema and Mitchell-Bruker, 2005;Wörman et al., 2006). These flow directions are, in turn, used to calculate topographic indices, such as the Topographic Wetness Index (TWI; Beven and Kirkby, 1979), that are widely used to describe the topographic controls on the spatial variation in wetness conditions. The topographically-based flow directions also allow the determination of the catchment boundaries based on the assumption that groundwater flow is predominantly parallel to the surface topography. Even though it is well known that the groundwater catchment boundaries may not be the same as the surface topography-based catchment boundaries, the surface topography-based drainage divide is usually used to determine the drainage area (i.e., catchment area) because data on the location of the groundwater divide is usually not available.
High-resolution elevation data also include many small-scale features that are unlikely to affect the direction of groundwater flow. In areas where the groundwater table is close to the surface, the flow pathways are highly affected by this microtopography (e.g., Bresciani et al., 2016b;Frei et al., 2010;O'Loughlin and Resources, 1986;van der Ploeg et al., 2012) but at better drained hillslope sites, the general flow direction at the hillslope or catchment scale is unlikely to be affected by the microtopography (Marklund and Wörman, 2011;Zijl, 1999). For the identification of the direction of slope parallel groundwater flow and the delineation of the catchment boundaries, preprocessing of highresolution DEMs has, therefore, become a common practice. It includes smoothing, aggregation or resampling of the DEM, as well as the removal of sinks (Leach et al., 2017;Lidberg et al., 2017;Lindsay and Creed, 2005;O'Neil et al., 2019). For example, it has been shown that smoothing improves the performance of a wetland classification model by effectively removing 'topographic noise' and small landscape features that would incorrectly be classified as wetlands (O'Neil et al., 2019). Aggregation (i.e., the use of a lower resolution DEM) can disturb the delineation of the stream network (e.g., Dehvari and Heck, 2013;Goulden et al., 2014;Lidberg et al., 2017;Mcmaster, 2002) but can also lead to more realistic groundwater flow directions (e.g., Gyasi-Agyei et al., 1995;Quinn et al., 1991;Walker and Willgoose, 1999;Wolock and Price, 1994;Zhang and Montgomery, 1994). Similarly, we expect that smoothing of small-scale topographic features that are unlikely to affect the groundwater table (either through actual smoothing or aggregation of the elevation data) can improve the simulation of the shape of the groundwater table and slope parallel groundwater flow directions.
Previous studies have looked at the effects of DEM resolution on the hydrological indices that are used for the characterization of surfaceand groundwater flow, but they mostly used grid resolutions from ~30 m (1 Arcsec) upwards (Le Coz et al., 2009;Stenta et al., 2017;Vieux, 1993). These resolutions are apt for large catchments but are often too Fig. 1. Map of the Krycklan catchment with the stream network, the 5 m contour lines and the location of outlets used in this study. The upper inset shows the location of the catchment in Sweden; the lower inset zooms into the focus area and shows also the 1 m contour lines (the focus area is located within the subcatchment that is called C6 in other studies (e.g. Laudon et al., 2013)).
coarse to capture the large variability in groundwater levels in small headwater catchments. Rinderer, van Meerveld and Seibert (2014) compared DEMs with a 2, 4, 6, 8, and 10 m resolution and considered the 6 m DEM most adequate because it reflected the main ridges and depressions but smoothed out the microtopography that they assumed had little effect on shallow groundwater flow directions. They, therefore, used this DEM to determine the topographic indices (e.g., TWI, curvature) that they related to the observed groundwater levels (Rinderer et al., 2014) and used to predict the groundwater level at unmonitored locations (Rinderer et al., 2019). Dehvari and Heck (2013) used 1, 2, 3, 4, 5, 8 and 10 m resolution DEMs to study the effect of DEM resolution on landscape attributes, such as ridges and flat areas, and concluded that the used coarser resolution strongly affected the extraction of streams. Similarly, Woodrow et al. (2016) used a 1, 5 and 10 m DEM and showed that contributing areas and stream lengths were significantly impacted by the choice of the DEM. However, no study has systematically investigated how the smoothing or aggregation of a DEM affects the estimated shallow groundwater flow directions or where in the landscape the calculated groundwater flow directions are most impacted by DEM smoothing or aggregation.
If the resolution of the DEM significantly affects the calculated groundwater flow directions, this also affects the location of the derived catchment boundaries. Previous studies have shown that catchment delineation based on almost unaltered (i.e., only minimal smoothing, aggregation, sink filling or channel burning for sink removal) topographic data can lead to either an overestimation or underestimation of the catchment area. The degree of error depends on the topography and geology of the catchment (Condon and Maxwell, 2015), as well as the resolution of the elevation data used for the delineation (Hinton et al., 1993;Woodrow et al., 2016).
In this study, we, therefore, focus on the effects of DEM-smoothing and -aggregation on the calculated directions of shallow groundwater flow in small headwater catchments and the derived catchment boundaries. Although groundwater flow is a three-dimensional phenomenon (Toth, 1963;Winter et al., 1999), and the deeper flow pathways of the nested groundwater system are important for the water balance and baseflow generation (Fan and Schaller, 2009;Oda et al., 2013;Welch et al., 2012), we focus here on the shallow groundwater flow pathways (typically within 6 m of the soil surface) in small (<1 km 2 ) headwater catchments. We assume that the groundwater table is a subdued copy of the topography and the main direction of shallow groundwater flow is parallel to the surface. We systematically smoothed high-resolution elevation data for the Krycklan catchment in northern Sweden to eliminate small-scale features from the DEM. This procedure corresponds, in principle, to reducing the resolution of the DEM. The idea is that after some degree of smoothing, the DEM of the surface topography better represents the groundwater table and flow directions (cf. Rinderer et al., 2014). If the DEM resolution has a significant impact on the calculated flow directions, the obvious follow-up question (which we do not address in this study) is whether there is an optimal resolution. However, as a first step, and in the absence of observed groundwater level data to determine the location of the groundwater table and flow directions, we focus here on quantifying the impact of DEM resolution (and thus information content) on the calculated groundwater flow directions and catchment areas for an area where shallow groundwater is the main source of runoff. More specifically, we address the following questions: • How much does the DEM resolution or DEM smoothing influence the calculated flow directions?
• Where in the landscape are the calculated flow directions most sensitive to the DEM resolution and information content?
• How much does the preprocessing of the DEM affect the estimation of catchment areas?

Study catchment
The 67.8 km 2 Krycklan catchment is located about 100 km northwest of Umeå in Sweden (see Fig. 1). It is located in the boreal zone and has a typical glaciated terrain. We selected the Krycklan catchment for this study because of the availability of a high-resolution (2 m) DEM and the relatively simple hydrogeological situation due to the deep till soils (Lindqvist et al., 1989). The groundwater tables are generally shallow (< 6 m from the surface and in most locations < 2 m). Subsurface flow pathways to the stream can be predicted based on the surface topography (Leach et al., 2017;Ploum et al., 2020). Shallow groundwater flow is the dominant source of streamflow (Lyon et al., 2012); in the riparian zone the groundwater level fluctuations are highly correlated with streamflow (Seibert et al., 2003).
The long-term annual average precipitation is about 614 mm/yr; the mean annual temperature is 1.8 • C (-9.5 • C in January, +14.7 • C in July). The mean annual runoff is 311 mm/yr, which means that about half of the precipitation goes to evapotranspiration (Laudon et al., 2013). The landscape consists mainly of forests, peatlands and lakes. The slopes in the catchment are not very steep, except for some incised stream banks (up to 10 m in the glacio-fluvial sediments). Human-made  (Lindsay, 2016a)). ditches can be up to 1 m deep and are primarily located in the headwater catchments with till soils. For more information on the Krycklan catchment, see Laudon et al., (2013).
For more detailed analyses and visualization of the results, we focused on a small area within the Krycklan catchment (named C6 in other studies in the Krycklan catchment). This focus area includes many characteristics and small-scale topographical features that are typical for the larger Krycklan catchment. This area is well known to the authors from field visits, which allowed for a better interpretation of the topographic features and the impact of smoothing or aggregation on the calculated flow directions.

General approach
The high-resolution DEM was modified in two different ways to obtain a lower resolution DEM: smoothing and aggregation. While the first approach assumes that the smoothed DEM might be a better representation of the groundwater table, the latter (aggregation) represents the situation where only coarser DEMs with a lower information content are available. We smoothed the DEM with two moving window filters: Gaussian and mean. Both methods employ a window of a defined pixel width that moves over the raster. The center pixel of the window is assigned a new value that is calculated either by weighing the values of the pixels in the window based on their distance according to a Gaussian curve with a defined standard deviation (Gaussian smoothing) or by taking the mean of all the pixels in the window (mean smoothing). We also aggregated the DEM (i.e., assigned the average value of all aggregated pixels to the new pixel) to obtain DEMs with a different information content. Aggregation can to a certain degree be considered a form of smoothing as it also leads to small scale features being averaged out. For the different DEMs, we determined the slope for each pixel. We assumed that this slope represents the direction of groundwater flow. We did this to determine where in the catchment the calculated flow directions are most sensitive to modifications of the topographic data and to determine how the catchment boundaries change when different DEMs are used.
Most steps of the workflow for the processing of the DEMs (Fig. 2) were implemented in Whitebox GAT (version 3.4 'Montreal', Lindsay, 2016a). The main steps and individual tools for these procedures are described below; additional references to the tools are provided in the supplementary material (Table S1).

Original DEM
The original 2-m DEM for the Krycklan catchment was based on LiDAR data obtained in August 2015. The LiDAR data has an average pulse density of 20 pulses/m 2 . The DEM was created by Gudrun Norstedt (Norstedt, 2017) and has been partially tested against ground observations during field visits (personal communication William Lidberg, Thomas Grabs, SLU). This DEM is referred to as the original DEM (orig) from here on.

Smoothed DEMs
For the smoothing of the original DEM, we used a Gaussian filter (Whitebox-Gaussian Filter, Lindsay, 2016a) and a moving window (arithmetic) mean filter with different window sizes (Whitebox -Mean Filter, Lindsay, 2016a). The window sizes were chosen to represent different degrees of smoothing of the surface (see example in Fig. 3), and therefore the assumed groundwater table. For the Gaussian filter, we  Fig. 1), showing the elevation for the original DEM with a 2 × 2 m resolution (orig, a) and examples of the smoothed and aggregated DEMs. Gaussian filter using a 3 × 3 pixels (gauss 3 , b) or 20 × 20 pixels (gauss 20 , e) standard deviation, Moving window filter using the mean value for a 3x3 pixel (mean 3 , c) or 21 × 21 pixels (mean 21 , f) window, and the aggregation of 3 × 3 (aggrg 3 , d) or 20 × 20 (aggrg 20 , g) pixels to obtain a lower resolution DEM. The field validated locations of the streams are shown in blue for reference. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) derived DEMs based on a weighted mean of the surrounding pixels, with the weights assigned by the Gaussian curve with a standard deviation of either 2 × 2 pixel, 3 × 3 pixel, 5 × 5 pixel, 10 × 10 pixel or 20 × 20 pixels. We refer to these DEMs as gauss n , where n refers to the number of pixels used for the calculation of the standard deviation that was used for the smoothing (2, 3, 5, 10 or 20). For the mean filter, we used the 3 × 3 pixel, 5 × 5 pixel, 9 × 9 pixel, 21 × 21 pixel window sizes because of the requirement for an odd number of pixels. These DEMs are named mean n , where, n again refers to the number of pixels used for the smoothing (3, 5, 9 and 21).

Aggregated DEM
We used aggregation to obtain lower resolution data (Whitebox -Aggregate, Lindsay, 2016a). For this, we aggregated the pixels (2 × 2 pixel, 3 × 3 pixel, 5 × 5 pixel, 10 × 10 pixel, 20 × 20 pixel) of the original DEM and assigned the mean value to all resultant pixels. Note that this aggregation differs from mean smoothing for which the mean is calculated for the center pixel of the window without changing the resolution of the raster (see Fig. 3). The aggregated DEMs are referred to as aggrg n , where n refers to the aggregated number of pixels (2, 3, 5 10 or 20). The number of pixels used for the aggregation was chosen to obtain data representing commonly available high-resolution DEMs (such as 5 m, 10 m, 25 m and ~30 m (1 Arcsec)). For example, Rinderer et al. (2017) used a 6 by 6 m DEM (which would represent a smoothing of 3 × 3 pixels if applied to our original DEM) to determine the flow directions and topographic indices that were correlated to the variation in groundwater levels in a pre-alpine catchment. Leach et al. (2017) resampled the DEM to 5 and 10 m (which would be a 2.5 × 2.5 and 5 × 5 pixel aggregation if applied to our original DEM resolution) to investigate where the groundwater table is close to the surface and the land is susceptible to damage by heavy machinery.
After calculating the flow direction and accumulated area (see below for more details), the aggregated rasters were resampled to match the extent and pixel size of the original and the smoothed DEMs. This means that the rasters resulting from the aggregation process nominally had the same resolution but a lower information content than the original DEM. In other words, the 'lower resolution' rasters allowed for pixel-by-pixel comparison with the other rasters, while de facto still having a lower resolution.

Hydrological conditioning of the DEMs
The original, smoothed and aggregated DEMs were hydrological conditioned prior to the calculation of the flow directions. Hydrological conditioning of a DEM means that the DEM is modified to remove any sinks, i.e., to ensure that for every pixel, there is a flow pathway to an outlet pixel with a strictly monotonically decreasing elevation. This conditioning is needed to ensure a realistic prediction of the locations of streams and channels , as well as wetland identification (O'Neil et al., 2019). The conditioning generally includes filling and or breaching of depressions (Lindsay, 2016a) to avoid problems with drainage area delineation due to small residual sinks and dams (Band, 1986;Lindsay, 2016b) that either naturally occur in the landscape or are created in the process of smoothing and aggregating the DEM. Lidberg et al. (2017) concluded that for the Krycklan area, a breaching algorithm leads to the most realistic results for stream and channel locations. We,  Figure S1. therefore, used a breaching algorithm (fast breaching algorithm in Whitebox (Lindsay, 2016a;Soille, 2004;Wang and Liu, 2006) for all the DEMs (i.e., original, smoothed, and aggregated DEM).

Flow direction analysis
For each of the 15 conditioned DEMs (one original, five Gaussian smoothed DEMs, four mean filter smoothed DEMs, and five aggregated DEMs), we calculated for each pixel the direction of the surface. We assumed that these surface gradients represent the shallow groundwater flow directions and, thus, refer to them as flow directions from hereon. We used three different methods to derive the flow direction for each pixel and each DEM: the slope aspect, the D8 (O' Callaghan and Mark, 1984), and D inf (Tarboton, 1997) algorithms. This resulted in 45 different flow directions for each pixel (15 DEMs multiplied by three different flow direction algorithms). Although algorithms that allow for multi-directional flow (e.g., Quinn et al., 1991;Seibert and McGlynn, 2007) are considered superior in certain situations, we decided not to use them here because of the difficulty in comparing the results when there are several directions per pixel.
The 45 flow directions for each pixel were compared using the circular variance (C V ) to determine the sensitivity of the calculated flow direction to the choice of DEM and flow direction algorithm. For this, we converted the flow directions from all 2x2 m DEMs (i.e., the original DEM and the outputs of the smoothing and aggregation steps) to radians. We then calculated the circular variance for each pixel using the 'raster' package (Hijmans, 2019) in R software (R Core Team, 2019). The circular variance C V is 1 minus the length of the mean resultant vector (ρ; see Eq. (1)) and indicates the clustering of the calculated flow directions, or a lack thereof.
The mean resultant length (ρ) was calculated using Eqs. (2)-(4) (after Pewsey, Neuhäuser and Ruxton, 2013), where x is the angle for the individual flow direction vectors i in radians and n is the total number of flow direction vectors (=45 in this case).
The value of the mean resultant length is one if all flow directions are similar and zero if they point in all directions (Kutil, 2012;Pewsey et al., 2013). For a more intuitive interpretation of the results (i.e., a larger value indicates more variation), we used the circular variance (C V ) (Eq. (1)) instead of the resultant mean length (Kutil, 2012). Small C V values indicate small variations in the flow directions and large values indicate a large variation (for examples, see Fig. 4 and Figure S1). For the interpretation of the results, we defined four classes of circular variance: low (C V < 0,026), middle (0.026 ≤ C V < 0.104), high (0.104 ≤ C V < 0.378) and very high (C V ≥ 0.378). These classes represent an approximately even distribution of flow directions within a range of 0 • − 45 • (low), 45 • − 90 • (middle), 90 • −180 • (high) and > 180 • (very high) (Fig. 4). These ranges were calculated based on a synthetic dataset of 45 evenly distributed flow directions over a range of 1 • to 360 • ( Figure S1).
A drawback of the resultant vector length and thus also the use of the circular variance (C V ) is that it does not provide clear information on the distribution of the flow directions. A value close to zero can, for example, indicate that the flow directions are evenly split in two opposite directions or that the flow directions are evenly spread into all directions (Pewsey et al., 2013). However, visual inspection of the calculated flow directions suggested that switches between two opposite directions were rare.
The variability in the flow directions (and thus the value of C V ) is affected by both the DEM and the flow direction algorithm. The differences in the calculated flow directions for a pixel for the three different flow direction algorithms were, however, almost always smaller than 45 • , and generally smaller than 20 • (i.e., within the low and middle variability category) ( Figure S2). The differences in the flow directions calculated for the three algorithms were, furthermore, larger for the original DEM than for the smoothed or aggregated DEMs, and, not surprisingly, differences in the flow directions calculated with the D8 and D inf algorithms were smaller than the differences in the flow directions calculated with either the D8 or D inf algorithm and the slope aspect algorithm (see Figure S2). In other words, the differences in flow directions due to the choice of the flow direction algorithm were small, and changes in the flow directions in the high and very high variability category are mainly caused by the choice of the DEM (i.e., due to the effects of smoothing or aggregation).
To determine the characteristics of the areas where the flow directions are highly sensitive to the choice of the DEM, we determined the Spearman rank correlation (r S ) between the circular variance and several topographic indices that were calculated for the original DEM: accumulated area, Topographic Wetness Index (TWI, Beven and Kirkby, 1979), distance to stream (downslope distance to stream; Lindsay, 2016a), downslope index (Hjerdt et al., 2004), and roughness index (Riley et al., 1999). Different landscape elements can have a similar topographic characteristic (e.g., both ridges and incised streambanks in the study area are steep). Therefore, we also determined the fraction of pixels with a particular combination of topographic index values with a high or very high circular variance. If a larger fraction of pixels with a certain combination of topographic index values (e.g., steep slope and large distance to stream for the ridges and steep slope but small distance to stream for the incised streambanks) has a high or very high circular variance, then the flow directions in this type of area are considered to be sensitive to the choice of the DEM.

Determination of the location of the drainage divides and accumulated areas
We also assessed the impacts of DEM smoothing and aggregation on the inferred position of the drainage divides. To assess the combined effect of potential shifts in the location of the drainage divides, we calculated the accumulated area (i.e., drainage area) for 40 points in the Krycklan catchment (see Fig. 1). For the selection of these 40 points, a ground-truthed stream network was combined with a map of the D8 based flow accumulation for the original DEM. The 40 points were chosen to represent the outlets of sub-catchments with accumulated areas of 10-12 ha, 120-140 ha, and 300-320 ha. We refer to these points as the outlets of the small, medium and large catchments, respectively.
We calculated the D8 flow accumulation for all DEMs and determined the pixel with the largest accumulated area within a 10 m radius from the original outlet position to consider potential shifts in the location of the simulated stream network for the smoothed and aggregated DEMs. More specifically, we used a snapping algorithm (based on Lindsay et al., 2008) to adjust the outlet position. To avoid that the snapping leads to a 'class jump', i.e., a point where the outlet no longer represents the original sub-catchment but is moved to a stream that drains a larger catchment, we chose only points for which there was no other stream with a similar or larger accumulated area within a 40 m radius for the original DEM. Even though the outlets were chosen based on the original DEM, the snapping process was applied for all DEMs (including the original DEM). For the comparison of the size of the calculated accumulated areas for the different DEMs for each outlet, we determined the relative area, which is the ratio of the accumulated area calculated for the smoothed or aggregated DEM and the accumulated area calculated for the original DEM.
The smoothing and aggregation process can result in unrealistic stream networks and catchment boundaries, particularly for small catchments (e.g., Goulden et al., 2014). With the snapping, we tried to account for small shifts in the stream network caused by the smoothing and aggregating (especially in cases where the stream position is on the border of two pixels). Still, there were cases where the smoothing and aggregation process caused the original stream to be shifted rather far away, or to no longer exist at all. This was particularly the case for the small catchments for which the stream was sometimes smoothed over, and the drainage area became part of a larger stream. We included these catchments in the analyses, even if their form and size clearly indicated that the smoothing or aggregation had resulted in an unrealistic catchment area. Excluding these cases would have been subjective in many cases because the changes were often gradual as the degree of smoothing/aggregation increased. Although all catchments were included in the statistics, we focused on the accumulated areas that were at least half but less than twice the size of the accumulated area derived for the original DEM (i.e., 0.50 ≤ normalized accumulated area < 2.0) because larger differences reflect larger structural changes in the stream channel network and a shift in the derived stream channel further away from the original outlet than the 10 m snapping distance. Such extreme changes in drainage area are more likely to be noticed by researchers working in these areas (although this is likely not the case when the algorithms are used for many catchments in a region) and could be filtered out by looking at the variation in the areas for each catchment for the differently smoothed DEMs (e.g., a very large variation could indicate a stream jump). Smaller changes due to the preprocessing of the DEM will be more difficult to detect but can still significantly influence water balance calculations.  Fig. 4 and Figure S1).

Variation in flow directions
For the majority of pixels, the calculated flow directions differed for the different DEMs and flow accumulation algorithms (Fig. 5). The flow directions tended to differ more from those calculated for the original DEM when the degree of smoothing or aggregation increased (see Fig. 6c and Fig. 7). The circular variance was low (C V ≤ 0.029) for only 28% of the whole Krycklan catchment (Fig. 5). For almost 24% of the pixels, the circular variance was high, and for about 6% of the pixels, the variance was very high (C V > 0.408). The calculated flow directions for the 0.2 km 2 focus area varied slightly more than for the entire 67 km 2 Krycklan catchment, but the distribution of C V was overall similar (Fig. 5).   Fig. 7. Boxplots showing the distribution of the absolute differences between the flow directions calculated for the adjusted DEMs and the original DEM for all different DEMs (Gaussian and mean smoothing and aggregation) and the three methods to calculate the flow directions (aspect, D8, Dinf). The box represents the 0.25 and 0.75 quantiles, the whiskers extend to the 0.1 and 0.9 quantiles, and the horizontal line within the box represents the mean. Figure S2 shows the difference in the calculated flow directions for the three methods used to calculate the flow directions for each DEM.

Table 1
The Spearman rank correlation (r s ) between the circular variance C V and the topographic indices (calculated for the original DEM) for the entire Krycklan catchment and the focus area. The calculated p-value were < 0.001 for all correlations due to the very large number of data points.  . 8. Proportion of the pixels with a certain distance to stream and local slope gradient (upper panels, a and b) or Topographic Wetness Index and local slope gradient (TWI, lower panels, c and d) for which the circular variance (C V ) was high or very high (C V > 0.104) for the entire Krycklan catchment (left-hand panels, a and c) and the focus, area (right-hand panels, b and d). The slope gradient, distance to stream and TWI were calculated for the original DEM. Figure S4 shows the number of pixels in each C V and slope gradient class for the entire Krycklan area.
The Spearman rank correlation (r S ) between slope (based on the original DEM) and the circular variance was −0.50 for the entire Krycklan catchment and −0.41 for the focus study area. The correlations between either the downslope index or the roughness index and the circular variance was also relatively high (Table 1). For most steep and very steep slopes (>40% based on the original DEM), the calculated flow directions were relatively stable (see Figure S4). For moderately steep slopes (up to 40%) and flatter areas, the probability of a high circular variance was larger. Still, the variation in flow directions could be either small or large (Fig. 8). For some moderate to steep slope areas with high TWI values (i.e., large accumulated area), however, the circular variance was very high (see Fig. 8c and d). These areas mainly represent the steeper incised streambanks, where the flow directions were impacted by a shift in the location of the stream channel due to the smoothing or aggregation. For the focus area where the streams are less incised, we did not find this effect (see Fig. 8a and b). The calculated flow directions for pixels near local water divides (low TWI and varying slope, Fig. 6a) also tended to vary more for the different DEMs (Fig. 8).

Catchment boundaries and accumulated areas
The smoothing and aggregation also affected the derived location of the drainage divides and, thus, the calculated accumulated areas for the selected outlet points (Fig. 9). More smoothing and aggregation generally led to smaller or larger accumulated areas (see Fig. 10 and Table S2). The accumulated areas of smaller catchments (as determined from the original DEM) were generally more impacted by the smoothing or aggregation than those of bigger ones (see Fig. 11). The most extremely smoothed DEMs tended to result in very small accumulated areas (Fig. 10) for the outlets of some of the small catchments. This was especially the case for Gaussian smoothing, but can also be seen for mean smoothing. In some cases, the shape of the accumulated area for an outlet was so different for the smoothed DEM that there was barely any overlap with the accumulated area calculated for the original DEM (Fig. 9a). Aggregation tended to lead to larger catchments ( Fig. 10) but there was no clear relation with catchment size (Fig. 11).

Effect of DEM smoothing and aggregation on the calculated flow directions
DEM smoothing or aggregation significantly affected the calculated flow directions for many pixels. High values of circular variance were mainly found in four different types of areas: i) areas where the local slope gradient differs from the main slope gradient, ii) flat areas, iii) ridges, and iv) the banks of incised streams. For areas where the local slope gradient differed from the general direction of the slope gradient, the flow directions shifted with increasing smoothing towards the "regional" trend ( Fig. 12a) because smoothing and aggregation make the topography more uniform and cause a flattening of these smaller local features. Examples of these types of areas are the esker and other small hills or bumps in the landscape (see #1 in Fig. 6a).
For flat areas, small changes in the elevation due to smoothing or aggregation can easily result in a change of the slope aspect (Fig. 12b). The circular variance values for these flat areas were usually medium to high (see, for example, #2 in Fig. 6a). This was particularly the case for larger flat areas in the valley bottom (high TWI and low slope gradient; Fig. 8c and d), but also for flatter areas near the ridges (large distance to stream and low slope gradient, Fig. 8a and b).
The circular variance values were generally also high for areas near the ridge, even for relatively steep slopes (low TWI or large distance to stream and steep slope; Fig. 8c). On the map, these ridges are often clearly visible as lines with high circular variance (e.g., #3 in Fig. 6a). The smoothing and aggregation of pixels at a ridge can cause a small shift in the location of the highest points (and thus the drainage divide) Fig. 9. The number of times that a pixel was considered to be part of the accumulated area for outlet points 1 (a) and 38 (b) for the different DEMs (the number of catchment layers is 15 in both cases but for outlet 1 there is no pixel where all layers overlapped which is why the scale only reaches up to 11). and a dramatic change in the calculated local flow direction (i.e., a shift from one side to the other; Fig. 12c). This effect was also observed for the incised streambanks (Fig. 8a, small distance to stream and steep slopes), where the smoothing could cause the stream to shift to a neighboring pixel (Fig. 12d).
We do not know which degree of smoothing or aggregation leads to the most suitable DEM to determine the groundwater table and groundwater flow directions. Field observations or modeling studies are needed to confirm which DEM, and thus which smoothed surface, most closely resembles the actual groundwater table. Still, this study highlights that smoothing of the DEM affects the calculated flow directions and the types of areas where the calculated flow directions are most affected by the smoothing. The results thus highlight the need to carefully consider the impacts of smoothing on the calculated flow directions and derived topographic indices, as well as the catchment area.
Most likely, there is not one single best DEM to describe the groundwater table. As the groundwater level drops, the groundwater table will not only be lower but also smoother (i.e., more subdued) than the surface topography (Bresciani et al., 2016a;Gleeson et al., 2011;Haitjema and Mitchell-Bruker, 2005;Marklund and Wörman, 2011;Winter, 2001). The systematic shift in the calculated flow directions for some areas, particularly areas where local features created a slope gradient away from the main hillslope gradient, could indicate that groundwater flow directions depend on the groundwater level. Several previous studies have shown that flow directions on hillslopes and in riparian zones can change dramatically as the groundwater level rises or falls (Rodhe and Seibert, 2011;van Meerveld et al., 2015;von Freyberg et al., 2014). If the degree of smoothing is a function of the depth to the water table, then at least for some locations in the catchment, different types of smoothing may be needed for different wetness conditions or seasons. Fig. 10. Cumulative frequency distributions of the normalized accumulated areas for the Gaussian smoothed DEMs (top), the mean moving window smoothed DEMs (middle), and the aggregated DEMs (bottom). Note that the axis is limited to areas that are at least half but less than twice the area calculated for the original DEM (i. e., 131 outliers in total are not shown). Also note the log scale for the normalized area.

Effect of DEM-smoothing and -aggregation on the location of drainage divides
The change in the calculated flow directions with smoothing or aggregation of the DEMs for the ridge sites causes a change in the location of the water divide and thus the surface topography-based catchment boundaries. Our results show that, particularly for small headwater catchments, the calculated drainage area depends strongly on how much the DEM is smoothed. For some DEMs and outlets, the drainage areas changed considerably (by much more than a factor of two) from those determined for the original DEM. This was most often the case for the outlets of the smaller headwater catchments but was, for some DEMs, also observed for the outlets of the medium and large catchments (Fig. 10). It was challenging to determine what caused these substantial changes in the drainage area. Still, visual inspection of the catchment boundaries and the DEMs suggested that in these cases, the DEM derived stream network changed considerably so that either a much larger area drained to the chosen outlet point, or more frequently there was no longer a stream that drained to the original outlet point because the stream channel had been smoothed out. This tendency for the streams of the smallest catchments to disappear when smoothing with the Gaussian or mean method should be taken into account, when determining catchment areas because it could lead to an underrepresentation of the catchment size of small headwater streams.
The effects of the aggregation on the relative catchment sizes are interesting because these results most closely emulate the impact of lower resolution elevation data (e.g., Gyasi-Agyei et al., 1995;Quinn et al., 1991;Walker and Willgoose, 1999;Wolock and Price, 1994;Zhang and Montgomery, 1994). For the more aggregated DEMs, there was a tendency for overestimation of the drainage areas, but this effect did not seem to depend on the (original) catchment size. Smoothing led to some very small "left-over" catchments as parts of the catchment were smoothed out and became part of other catchments. Aggregation did not lead to (as many of) these "left-over" catchments because of the larger cell-size. However, in some cases, it led to the merging of catchments and thus the frequency of large catchment areas increased.
The differences in the derived catchment areas are important, especially when considering how they will affect the water balance for these catchments. While field knowledge will help spot a very incorrect location of the drainage divide, small changes in its location and the derived catchment areas will be more difficult to detect and correct. The results from this study suggest that without any field confirmation, the computed catchment areas for small catchments have to be considered as very uncertain. A plausibility check is not possible when DEMs are Fig. 11. The normalized accumulated area for each outlet for the different DEMs (dots, jittered for better clarity) and corresponding boxplots. Note that the areas that are less than half and more than double the area calculated for the original DEM are not shown and not included in the boxplots, but are indicated using grey arrows, with the numbers displaying the number of outliers (131 outliers in total, i.e., 20 % of all values). The outlets were sorted according to the accumulated area calculated for the original DEM (smallest at the bottom). Note the log scale for the normalized area. The box represents the 0.25 and 0.75 quantiles, the whiskers extend to the 0.10 and 0.90 quantiles, and the line within the box represents the mean. used to automatically delineate many headwater catchments. Water balance data are usually not available for these cases either so that it will be difficult to detect the incorrect catchment areas. Potentially, the methodology used here (i.e., calculating the drainage area for different DEMs) could provide information on the uncertainty of the calculated catchment areas.
A shift in the location of specific landscape features, such as ridges, with increasing smoothing of the surface, could also be an indication of groundwater subsidies to or from other catchments (Tóth, 1962). Indeed systematic water balance differences have been observed between adjacent catchments in this landscape (Karlsen et al., 2016), as well as elsewhere (e.g., Oda et al., 2013). Our methodology can be helpful to identify areas of particular interest when deciding on the placement of wells for calculating groundwater gradients, which can be used to gain more insights into runoff generation or to determine the boundary of the groundwater divide.

Applicability for other regions
Our study focused on the groundwater table and assumed that groundwater flow is predominantly parallel to this surface (or horizontal, cf. the Dupuit-Forcheimer assumption). We, furthermore, assumed that the glacial till soils are laterally homogeneous in terms of hydraulic properties and the groundwater system is topography driven (Gleeson et al., 2011) so that the smoothed surface topography can be used to approximate the shape of the groundwater table. More sophisticated methods for the calculation of the groundwater table exist (e.g., Wörman et al., 2006;Zijl, 1999), but they are also often much more complex and require more in-depth knowledge of the investigated system. Although the shallow groundwater flow pathways are most important for runoff generation during rainfall or snow melt events in small catchments (e.g., Lyon et al., 2012), deeper (non-slope parallel flow) occurs as well and can be important for the generation of baseflow (Tóth, 1962;Winter et al., 1999). Even at Krycklan, the groundwater system is a nested system with shallow and deeper groundwater flow pathways (Kolbe et al., 2020). These deeper pathways can also cause inter-catchment groundwater flow (Fan and Schaller, 2009). The use of surface topography based estimates of groundwater flow directions and the drainage area is problematic if inter-catchment groundwater flow is important (Fan, 2019;Käser and Hunkeler, 2016;Oda et al., 2013;Welch et al., 2012). In other cases, only a part of the surface catchment may contribute to streamflow (e.g., Godwin and Martin, 1975;Liu et al., 2020;Shook et al., 2021).
The focus in this study was on small-scale boreal headwater catchments with shallow groundwater in relatively homogeneous till soils, but the results of this study are assumed to be applicable for other landscapes with relatively uniform and deep unconsolidated material and a water table that is relatively close to the surface. We assume that our approach and results are also applicable for landscapes with a more pronounced (i.e., steeper) topography, as long as the unconfined aquifer is still relatively homogeneous. In mountainous terrain and larger scales, the regional flow pathways may be more significant (e.g., Gleeson and Manning, 2008;Welch and Allen, 2012) so that the assumption of flow parallel to the surface is no longer valid. However, even in mountainous terrain, the near-surface flow component generally dominates unless the anisotropy is very large (Gleeson and Manning, 2008). Direct extrapolation to very flat areas may be problematic because the smoothing of local topographic features, such as the esker (Fig. 12a), is less likely where these features are not present but the sensitivity of flow directions to changes in the topographic information (Fig. 12b) is likely similar. However, for flat areas and wetlands, where the groundwater table is very close to the surface, groundwater flow may be driven by the microtopography (Frei et al., 2010;van der Ploeg et al., 2012), so that the smoothing would not lead to a better representation of the groundwater table and flow directions.
The results are most likely not applicable for areas where the depth of the soil or unconsolidated material is more variable, or where bedrock, fractures or impermeable layers have a large effect on the groundwater flow directions. Where the soils are shallow but a lower permeability layer is located close to the surface, the topography of this layer may determine the direction of groundwater flow. For example, van Meerveld et al. (2015) showed that in the Panola catchment in Georgia, U.S., the groundwater flow directions were more similar to the direction of the bedrock surface than the ground surface when the water levels were low, but more closely resembled the ground surface when the water levels were high.

Conclusions
In humid climates, the groundwater table is often assumed to be a replicate of the surface topography. This allows the use of digital elevation data to estimate the shape of the water table and the groundwater Fig. 12. Schematic drawings of the four topographic features for which the flow directions differed significantly for the different DEMs, leading to a high circular variance : a) areas where the local aspect differs from the general (i.e., regional) slope, b) a flat area, c) a ridge, and d) an incised stream. Original DEM (solid line) on top, smoothed DEM (dashed line) below. The arrows indicate the flow directions (solid for original DEM and dashed for the smoothed DEM) for the point of interest, indicated by the light grey vertical line. flow directions. However, some smoothing of high-resolution elevation data is needed to represent the groundwater table, as small-scale topographic features probably do not affect the groundwater flow. To investigate how the degree of smoothing of digital elevation data affects the calculated flow directions and where in the catchment the calculated flow directions are most affected by the chosen degree of smoothing, we generated different DEMs, representing different degrees of smoothing and aggregation. For each pixel, we then determined how the degree of smoothing affects the direction of the surface and, therefore, the assumed groundwater flow direction.
For almost a third of the study catchment, the calculated flow directions varied considerably for the different DEMs and the circular variance was either high or very high. For >70% of the catchment, the circular variance was at least as high as in a situation where the flow directions are uniformly spread over a 45 • range (medium circular variance). The calculated flow directions were most affected by smoothing in areas where the slope gradient was <40%. Smoothing or aggregation caused the flow directions to change particularly much for flatter areas with a high TWI and areas with local slopes in the opposite direction of the dominant slope gradient. These changes or uncertainties in the calculated groundwater flow directions also impact the derived locations of the drainage divides and the calculated accumulated areas. The effect was particularly larger for small headwater catchments (<12 ha), for which the smoothing often led to an underestimation of the catchment area.
While the concept of the groundwater table in humid climates being a subdued copy of the surface topography is well established, the question remains how subdued this copy is and which degree of smoothing or aggregation best reflects the actual groundwater table. Two approaches could be used to address this question: groundwater flow modelling or field observations of groundwater levels (and thus gradients). Our results demonstrate the need for such studies to determine the degree of smoothing of the DEM that yields the best approximation of the groundwater table.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.