Impacts of grid resolution on surface energy fluxes simulated with an integrated surface-groundwater flow model

The hydrological component of the Terrestrial Systems Modeling Platform (TerrSysMP), which includes integrated surface-groundwater flow, was used to investigate the grid resolution dependence of simulated soil moisture, soil temperature, and surface energy fluxes over a subcatchment of the Rur, Germany. The investigation was motivated by the recent developments of new earth system models, which include 3-D physically based groundwater models for the coupling of land–atmosphere interaction and subsurface hydrodynamics. Our findings suggest that for grid resolutions between 100 and 1000m, the non-local controls of soil moisture are highly grid resolution dependent. Local vegetation, however, strongly modulates the scaling behavior, especially for surface fluxes and soil temperature, which depends on the radiative transfer property of the canopy. This study also shows that for grid resolutions above a few 100m, the variation of spatial and temporal patterns of sensible and latent heat fluxes may significantly affect the resulting atmospheric mesoscale circulation and boundary layer evolution in coupled runs.


Introduction
In recent years, a growing number of earth system modeling platforms attempted to include physically based hydrological models, with lateral flow and groundwater surface water interactions, to study the linkages between landatmosphere and subsurface hydrodynamics (e.g., Anyah et al., 2008;Maxwell et al., 2011;Shrestha et al., 2014;Butts et al., 2014;Larsen et al., 2014). These studies show that the inclusion of groundwater dynamics improves the simulated spatial variability in root zone soil moisture and groundwater table depth, and shows the potential for improved forecasts of the whole terrestrial system. However, as soon as one moves from column-based land surface models to 3-D models with lateral flows, a new dimension of spatial complexity is added where scaling issues become highly relevant (Becker and Braun, 1999). This is mainly due to the introduction of non-local controls on soil moisture patterns (e.g., patterns of soil moisture dominated by lateral fluxes of surface and subsurface flow) as earlier identified by Grayson et al. (1997), which also depend on grid resolution. For spatial extents of 100 km and above atmospheric models are still run at grid resolutions ≥ 1 km due to computational limitations or physical parameterizations, and hydrological models coupled the atmospheric models are usually run at similar grid resolutions, which may however be inadequate to correctly simulate subsurface flow. Hyper-resolution models have already been suggested by, e.g., Wood et al. (2011), while Beven and Cloke (2012) have suggested the need for spatial-scale-dependent subgrid-scale parameterizations to adequately simulate soil moisture variability.
In reality, catchments exhibit variability and heterogeneity at a range of scales (Blöschl and Sivapalan, 1995), while in numerical models, the variability of soil moisture, soil temperature and surface fluxes can only be controlled by the heterogeneity at the chosen grid resolution. Previous studies with offline hydrological models have shown that the aggregation of topography to coarser grid resolution (e.g., 1 km) has a strong impact on the water balance (e.g., Zhang and Montgomery, 1994;Kuo et al., 1999;Vivoni et al., 2005; 4318 P. Shrestha et al.: Impacts of grid resolution on surface energy fluxes Bormann, 2006;Herbst et al., 2006;Giertz et al., 2006;Dixon and Earls, 2009;Sciuto and Diekkrueger, 2010;Sulis et al., 2011). Many of these studies focused primarily on the spatial-scale-dependent behavior of catchment discharge, groundwater table depth and catchment mean soil moisture. A better understanding of the spatial-scale dependence of the simulated patterns of soil moisture, temperature and surface fluxes is required, however, when such models are coupled to atmosphere models. In an idealized setup using a mosaic approach, Shrestha et al. (2014) demonstrated the importance of subgrid-scale topography on topographically driven surface-subsurface flow for land-atmosphere interactions and stressed the importance of accurate simulation of spatiotemporal variability of surface fluxes for the evolution of the terrestrial system as a whole.
The aim of this study is to examine the effects of resolution-dependent model heterogeneity using the Terrestrial System Modeling Platform (TerrSysMP, Shrestha et al., 2014;Gasper et al., 2014;Sulis et al., 2015) on the variability of modeled soil moisture, soil temperature and surface fluxes in a temperate climate when no subgrid-scale parameterizations are used. The rest of the manuscript is organized as follows: Sect. 2 describes the modeling tool used for the study; experiment design and setup of the catchment are discussed in Sect. 3. Topography heterogeneity analysis is presented in Sect. 4, while results and discussions are presented in Sect. 5, and conclusions in Sect. 6.

Modeling tool
The hydrological component of TerrSysMP consists of the NCAR Community Land Model CLM3.5 (Oleson et al., 2008) and the 3-D variably saturated groundwater and surface water flow code ParFlow (Ashby and Falgout, 1996;Jones and Woodward, 2001;Kollet and Maxwell, 2006;Maxwell, 2013). The two models ( Fig. 1) are coupled using the OASIS3 external coupler (Valcke, 2013). In the sequential information exchange procedure, ParFlow sends the updated relative saturation (S w ) and pressure ( ) for the top 10 layers to CLM. In turn, CLM sends the depthdifferentiated source and sink terms for soil moisture (top soil moisture flux (q rain ), soil evapotranspiration (q e )) for the top 10 soil layers to ParFlow (see Fig. 1). A more detailed description of the coupling can be found in Shrestha et al. (2014). In this study, the hydrological component of TerrSysMP is decoupled from its atmospheric component and forced with spatially distributed atmospheric forcing data at 2.8 km spatial resolution and hourly temporal resolution (air temperature (T ), wind speed (U ), specific humidity (QV), total precipitation (Rain), pressure (P ), and incoming shortwave (SW) and longwave (LWdn)) from COSMO-DE (Baldauf et al., 2011) analysis data of the German Weather Service (DWD). . OASIS3 is the driver of the component models for land surface (CLM) and subsurface (ParFlow). The configuration file for OASIS3 prescribes the end-point data exchange between the component models in a sequential manner. The variables exchanged between the two models are relative saturation (S w ), soil pressure head ( ), top soil moisture flux (q rain ), soil evapotranspiration (q e ), air temperature (T ), wind speed (U ), specific humidity (QV), total precipitation (Rain), pressure (P ), and incoming shortwave (SW) and longwave (LWdn).

Numerical experiment design
The model was set up for a sub-catchment of the Rur River (TR32 test bed site; Vereecken et al., 2010;Simmer et al., 2014), on the northern foothills of the lower Eifel mountain range with an approximated drainage area of 325 km 2 (Fig. 2a). The sub-catchment encompasses the tributary Wehebach, which merges with the River Inde. The elevation in the model domain reaches from 50 to 600 m from north to south (Fig. 2b), with mostly agricultural crops (c1n) near the foothills and needleleaf evergreen trees (nle) and broadleaf deciduous trees (bld) along the sloping terrain (Fig. 2c). The urban areas are represented in the model as agricultural canopy (c1f) with a fixed leaf area index (LAI = 0.6). Topography and land use are based on the 90 m resolution Shuttle Radar Topography Mission (SRTM) data and 15 m resolution data available from the TR32 database (Waldhoff, 2012), respectively.
The SRTM data were aggregated to 120, 240, 480 and 960 m horizontal grid resolution for the model domain setup (Table 1) by first interpolating the 90 m topography to 120 m using bilinear interpolation before aggregating to the coarser resolutions. Land use was aggregated by specifying only the dominant plant functional type (PFT) at the coarser resolution. The chosen grid resolutions are within the limits where a positive spatial autocorrelation of the topographic index exists (Cai and Wang, 2006), and roughly cover the range of grid resolutions between large-eddy simulation (LES) and mesoscale atmospheric modeling. For fully coupled mesoscale modeling with grid resolutions ≥ 1000 m (atmosphere component model), the above selected grid resolutions can be used for the hydrological component in TerrSysMP with the mosaic approach to better resolve the heterogeneity of topography, land use and geology.
The model setup for all grid resolutions used the same 10 vertically stretched layers (2-100 cm from top to bottom) followed by 20 constant depth levels (135 cm) extending to 30 m below the land surface. A uniform soil texture was used for this study by keeping the soil parameters spatially constant. The subsurface parameters were set as follows: saturated hydraulic conductivity, K s = 0.00034 m h −1 ; van Genuchten parameters, α = 2.1 and n = 2.0 m −1 ; and porosity, φ = 0.4449. This removes any impact of soil heterogeneity on non-local controls of simulated soil moisture variability and scaling of soil hydraulic properties. The soil moisture profile for all setups was initialized with a horizontally homogeneous hydrostatic pressure head and a water table depth at 5 m from the surface. The soil temperature was also initialized horizontally homogeneous with a uniform temperature of 10 • C for all levels. A time step of 3600 s was used, and the simulation was integrated for 6 years using hourly atmospheric forcing from COSMO-DE analysis data. The same atmospheric forcing data for the year 2009 were used recursively for the 6 years. The model outputs were averaged over 5 days for the analysis. The NCAR Command Language (NCL; NCAR, 2013) was used for data analysis.

Topography heterogeneity analysis
Heterogeneity analyses of the topography at different grid resolutions are summarized in Table 2. The profile and plan curvature represents the flow acceleration and convergent/divergent flow, respectively. The profile curvature is parallel to the direction of maximum slope and a positive/negative value indicates that the surface is upwardly concave/convex at the grid cell. The plan curvature is perpendicular to the direction of the flow and the positive/negative value indicates that the surface is convex/concave to the side at the grid cell (see https://resources.arcgis.com for pictorial descriptions). The distributions of plan and profile curvature change with coarsening of grid resolution. The plan curvature is negatively skewed at small grid resolution, and the skewness decreases with the coarsening, while skewness changes for the profile curvature are negligible. However, the kurtosis of both the plan and profile curvatures decreases exponentially with higher exponential power for plan curvature. This is also qualitatively visible in the streamline maps of D4 flow direction and the local slopes for the sub-catchment at the different grid resolutions (Fig. 3a, b, c, d): the aggregation of topography results in smoothing of slope magnitudes and the filtering of small-scale convergence and divergence zones. Without sub-grid-scale parameterization, this spatial filtering will impact lateral flow and simulated mean grid cell soil moisture distributions . Similar to the findings of Quinn et al. (1995), grid coarsening also affects the location of the water divides and makes it difficult to accurately delineate the catchment contributing area especially for d480 and d960.

Results and discussions
The simulated unsaturated storage (S unsat ) at different grid resolutions for the sub-catchment showed different temporal evolutions (Fig. 4). The unsaturated storage was normalized by the modeled sub-catchment area to account for differences in catchment size at different grid resolutions. The increase in S unsat during the first year reflects the adjustment of the subsurface storage from the horizontally homogeneous hydrostatic initial condition, with a ground to water table depth of 5 m. In the first half of this year, the finer grid resolutions adjust faster due to a more efficient drainage mechanism. In the second year, S unsat starts to decrease gradually, reaching a quasi-equilibrium in all simulations in the fifth year. This steady-state value range is lower for the coarser grid resolutions, caused by a higher groundwater table, which results from a less efficient drainage combined with higher infiltration at the lower resolutions. This is also illustrated in Fig. 4b, which shows that the average annual unsaturated storage (S unsat t ) in the sixth year is concurrent with the average slope of the sub-catchment. The decrease in average catchment slope with resolution coarsening is, however, also accompanied by decreasing plan and profile curvature kurtosis (see Table 2). These results are consistent with, e.g., Kuo et al. (1999), who related the increase in average soil moisture contents with grid coarsening to decreasing slope gradient and curvature variations. Sulis et al. (2011) also explained catchment wetness in terms of storage and ground to water table depth via decreasing local slopes and plan curvature variations due to aggregation effects. Different S unsat t with grid resolutions reflects different spatial soil moisture variability, which in turn influences simulated land-atmosphere interactions. Figure 5 shows the distribution of average top 10 cm relative soil moisture (S w ), average top 10 cm soil temperature (T soil ), sensible heat flux (SH) and latent heat flux (LH) for time periods when the model exhibits strong coupling with the atmospheric forcing. We assume strong coupling when the 5-day mean incoming solar radiation exceeds 128 Wm −2 , which corresponds to the period from April to September. We further filter for different PFTs to analyze the linkages between local vegetation and the non-local controls of soil moisture patterns with grid coarsening. Their respective temporal and catchment-averaged values (S w x,t , T soil x,t , SH x,t and LH x,t ) are indicated as solid markers (Fig. 5) and summarized in Table 3 for the discussion below. For relative soil moisture, grid coarsening leads to a sharp decrease in interquartile range and an increase in the 25 % quartile, median and mean value, i.e., reduced variability and higher mean simulated soil moisture. This is true for all PFTs. While no significant loss in interquartile range with reduced resolution is observed for temperature and the turbulent fluxes, a clear PFT-dependent scaling with grid coarsening exists, especially for crop PFTs and the crop PFT with a fixed low LAI. Average relative soil moisture increases by 30 % for trees and 23 % for crops when coarsening the grid resolution from 120 to 960 m. The difference between the different PFTs is partly related to their spatial location: trees are mostly located in steeper terrain compared to crops (according to the distribution of slopes for the different PFTs, not shown here). For trees, grid coarsening leads to lower T soil x,t by 0.6 and 0.3 • C for needleleaf evergreen trees and broadleaf deciduous trees, respectively, while for crops, T soil x,t is lowered by almost 1 • C. Thus forested grid cells exhibit higher grid resolution sensitivity for soil moisture and lower sensitivity for soil temperature, while grid cells with crops show the inverse. These findings can be attributed to the PFT-specific transmissivity for solar radiation, or the partitioning of absorbed solar radiation by vegetation and the ground. Needleleaf evergreen trees (nle) absorb more and transmit less solar radiation to the ground compared to broadleaf deciduous trees (bld), while the crop type with variable LAI (c1n) absorbs more solar radiation and transmits less compared to the constant LAI type (c1f). The PFT-specific optical parameters, including albedo and amplitudes of seasonal LAI, control the variability in solar radiation partitioning. This directly modulates the partitioning of sensible and latent heat flux and constitutes a local vegetation control. For all grid resolutions, the inter-PFT flux differences are much higher than the intra-PFT differences due to different grid resolutions. In general, grid cells with crops are more sensitive to grid resolution changes. A 60 % decrease in SH x,t and a 35 % increase in LH x,t are observed for c1f, with the grid coarsening from 120 to 960 m. SH x,t decreases by 30 % and LH x,t increases by 16 % for c1n. For forest, SH x,t decreases only by about 20 and 11 %, while LH x,t increases only by 8 and 11 % for bld and nle, respectively. In general, the latent heat flux increases and sensible heat flux decreases for all PFTs for coarser grid resolutions, while amplitudes of change are PFT-dependent, which suggests a local control by the PFT. This also suggests that the non-local control of soil moisture, which is affected by change in grid resolution, also controls the partitioning of surface energy fluxes, especially for PFTs transmitting more solar radiation to the ground. The average relative soil moisture for the subcatchment in this study was relatively wet; for drier regimes, the effects of model grid resolution on surface fluxes may be stronger. The above findings bear important consequences in terms of spatial variability of surface fluxes, which may affect the evolution of the atmospheric boundary layer, induce mesoscale circulations, and even lead to the formation of clouds and precipitation (Avissar and Schmidt, 1998;Baidya Roy and Avissar, 2002). This spatial variability is illustrated using a cross section along the sub-catchment. Figure 6a shows the spatial heterogeneity of topography and land use along cross-section AA indicated in Fig. 2 for different grid resolutions. The difference in variability of topography is not visible along the cross section due to the plot scale, but the effect of grid resolution on the simulated average top 10 cm soil moisture S w t is obvious in Fig. 6b. With grid coarsening, the finer scale plan and profile curvatures are filtered, which reduces drainage efficiency, increases the simulated grid cell soil moisture, and damps the spatial variability. This effect is more pronounced from d240 to d480 and d960 than from d120 to d240. Quantitatively, the spatial average and standard deviation for S w t along the cross-section AA are 0.70±0.16, 0.75±0.17, 0.90±0.10 and 0.93±0.04 for d120, d240, d480 and d960, respectively, which clearly shows the increase in soil wetness and decrease in variability with coarsening with largest changes between d240 and d480. The strong scaling behavior of soil moisture also modulates the partitioning of surface energy fluxes. Figure 6c shows the annual average Bowen ratio (ratio of sensible to latent heat flux) along cross-section AA ; its profile matches well with the PFT profile, indicating local vegetation control on the spatial pattern of surface flux partitioning. This could be partly enhanced due to the wet soil condition in the sub-catchment for the simulated period. However, the change in the non-local control of soil moisture with grid resolution also contributes to the profile of surface flux partitioning visible as a perturbation in the amplitudes of the Bowen ratio for d960. Some perturbations are also caused by different PFTs for finer grid resolution. Table 4 summarizes the statistics of the Bowen ratio profile along the cross section and shows Table 4. Annual average Bowen ratio and standard deviation along the cross-section AA for the land-use classes nle (needleleaf evergreen tree), bld (broadleaf deciduous tree), c1n (crops with seasonal LAI), and c1f (crops with fixed LAI).
PFT 120 m 240 m 480 m 960 m nle 1.56 ± 0.16 1.35 ± 0.22 1.22 ± 0.12 1.24 ± 0.07 bld 0.82 ± 0.08 0.72 ± 0.08 0.60 ± 0.09 0.59 ± 0.08 c1n 0.46 ± 0.13 0.50 ± 0.15 0.36 ± 0.13 0.28 ± 0.07 c1f 0.67 ± 0.28 0.53 ± 0.38 0.17 ± 0.11 0.18 ± 0.13 that its dominant pattern is strongly controlled by the PFT pattern. For trees, nle has a higher Bowen ratio (> 1) than bld (< 1), which is mainly due to differences in plant physiological properties, consistent with observations (Baldocchi et al., 1997). For crops, c1f has a higher Bowen ratio compared to c1n, which is due to the LAI difference. For both trees and crops, the Bowen ratio in general decreases with grid coarsening. Coarsening from d120 to d960 decreases the Bowen ratio by 20, 28, 39, and 73 % for nle, bld, c1n and c1f, respectively. The most significant change is found for crops with low LAI. Radiation absorbed by the ground plays a significant role in amplifying/attenuating the grid resolution dependence of surface flux partitioning. Again, it has to be mentioned that this statement may be valid only for wet regimes. Figure 7 shows the scatterplot between the Bowen ratio and average top 10 cm relative soil moisture along crosssection AA over the averaged time period filtered for different PFTs to illustrate the PFT-related dependence of the Bowen ratio on relative soil moisture. Large scatter for d480 and d960 m is found for S w < 0.7, while for d120 and d240, large scatter is only found for S w ≤ 0.4. Thus, the Bowen ratio distribution shifts with grid resolution. A linear regression between Bowen ratio and relative soil moisture gives a firstorder estimate of the Bowen ratio dependence on relative soil moisture. Tree canopies exhibit more variability in the Bowen ratio than crops. For trees, nle exhibits stronger scaling behavior with relative soil moisture than bld. For crops, c1f exhibits stronger scaling behavior with relative soil moisture than c1n. The results again show that crops with low LAI have a stronger influence on flux partitioning with grid resolution.
We also evaluated the grid resolution effects on fluxes using the mosaic approach by aggregating the simulated surface sensible and latent heat fluxes for d120, d240, and d480 to 960 m resolution. Figure 8a shows the time series of the 5day average of sensible heat flux along cross-section AA for d120. It shows the strong seasonal cycle of sensible heat flux that correlates with the seasonal cycle of net radiation (not shown here), and also the strong gradient along AA , owing to the different canopy cover. The differential heating along AA potentially generates mesoscale boundary layer circulations embedded in the local topographic circulation, whose strength would also depend on the mean heating rates of the catchment and the synoptic wind strength as indicated in many previous studies (Lemone et al., 2002;Baidya Roy and Avissar, 2002;Grossman et al., 2005). Figure 8b-d shows the difference in the time series of SH along cross-section AA between the coarser grid resolutions and d120. With grid coarsening the amplitude of the SH difference increases, suggesting an overall decrease in SH as also observed from the bulk quantities. Some d960 grid cells also exhibit an increase in SH mainly due to the changing dominant PFT, when the subgrid cells consist of trees and crops. Similarly, Fig. 9a shows the time series of 5-day average latent heat flux along cross-section AA for d120, which also exhibits a strong seasonal cycle and a strong gradient along AA . The difference in time series of LH along AA for coarser resolutions with respect to d120 shows the sharp increase in latent heat flux. The decrease and increase in SH and LH, respectively, are particularly high between 480 and 960 m resolution. Thus, when using coarser grid resolutions for the hydrological model, coupled to the atmospheric model, the simulated boundary layer would be relatively moister and cooler.

Summary and conclusions
This study was motivated by recent efforts in including physically based hydrological models in earth system models both for seasonal-scale and climate studies that would allow  for examination of the linkages between land-atmosphere and subsurface hydrodynamics. The hydrological component of the newly developed TerrSysMP was used over a sub-catchment of the Rur at grid resolutions encompassing roughly the spatial scales between LES and mesoscale atmospheric models to quantify the effect of grid resolution on simulated soil moisture, soil temperature and surface fluxes.
The terrain analysis of the sub-catchment showed the expected smoothing of slopes and filtering of the profile and plan curvatures with grid coarsening. This grid resolution has a strong effect on the non-local controls of soil moisture simulated by the model, while the local vegetation exerts a strong modulation on the transfer of the grid-resolutiondependent soil moisture variability to soil temperature and surface fluxes. In this study, soil moisture beneath forests was found to decrease more than beneath crops due to the location of forests over steeper slopes. However, due to the plant physiological properties affecting the transmissivity of solar radiation, crops lead to a higher grid resolution dependence than trees in terms of soil temperature and surface fluxes. For crops, the magnitude of the LAI was also found to have a strong effect on the scaling behavior of surface fluxes. This non-linear scaling behavior of the energy balance with respect to grid resolution can alter the spatial and temporal pattern of simulated surface fluxes. Larger differences were especially observed when moving from d480 to d960. These dependencies can induce or weaken mesoscale circulations and the ensuing boundary layer evolution when using coupled simulations. Using an idealized setup, such atmospheric feedback effects for a rainfall-runoff process (runoff due to excess infiltration, creating wet and dry patches), with varying grid resolution, were shown earlier by Shrestha et al. (2014). Here the heterogeneity at different resolutions was introduced in terms of topographic slopes while keeping a homogeneous land cover (crop). In this study, the inclusion of convergent zones at finer grid resolutions compared to coarser grid resolutions, where they are usually filtered out by topography smoothing, enhanced the overland flow, thereby reducing the infiltration and the mean soil moisture content while extending the downslope moist area. This extended downslope moist area increased the extent of the moist patch compared to the coarser run, thereby lowering the Bowen ratio along the extended patch, and also increasing the extent of the downdraft region in the atmospheric boundary layer. Thus, the dependence of the simulated surface-subsurface physical processes on grid resolution potentially also affects the local circulation and atmospheric boundary layer evolution. However, fully coupled real data simulations including the atmosphere at different grid resolutions are required to further improve our understanding of the land-atmosphere feedbacks.
The study was limited to grid resolutions from 120 to 960 m. One could argue that even the 120 m resolution is not sufficient for hydrological models, and much finer resolutions (≤ 30 m) are needed (e.g., Kuo et al., 1999). We acknowledge the limitation in this study, but finer resolutions challenge currently available computation resources and also the convergence of 3-D integrated surface groundwater mod-els. We realize that sub-grid-scale parameterizations along with resolutions of approximately 100-200 m would be sufficient for coupled simulations. The interaction between surface and groundwater in ParFlow is simulated using a 2-D shallow overland flow equation as an upper boundary condition (Kollet and Maxwell, 2006), instead of the commonly applied conductance concept. Particularly for the surfacegroundwater interaction, this upper boundary condition and the Darcy flux in the horizontal direction are affected by the spatial filtering of the terrain curvature in the coarsening of the grid resolution. This results in an increase in infiltration, a reduction of lateral flow and a shallower groundwater table. So, there is a need for a scale-dependent subgrid parameterization for the upper boundary condition and for the lateral flow in the groundwater model used in this study. For the upper boundary condition, e.g., the concept of a fractional saturated area (parameterized by the subgrid-scale distributions of topographic index and groundwater table depth) is widely used to control surface runoff and hence infiltration in many 1-D land surface models (e.g., Niu et al., 2005). For the subsurface, e.g., Niedda (2004) proposed an amplification/upscaling of hydraulic conductivity to compensate for the reduction in the hydraulic gradients using the information content of terrain curvature to produce similar lateral flow. In general, the topography heterogeneity analysis and the soil moisture data presented in this study do provide the necessary data to investigate such parameterizations for the upper boundary condition and lateral subsurface flow. Future studies will involve developing such robust scale-dependent subgrid parameterization for the 3-D physically based groundwater model.