Presentation and discussion of the high-resolution atmosphere–land-surface–subsurface simulation dataset of the simulated Neckar catchment for the period 2007–2015

Coupled numerical models, which simulate water and energy fluxes in the subsurface–land-surface– atmosphere system in a physically consistent way, are a prerequisite for the analysis and a better understanding of heat and matter exchange fluxes at compartmental boundaries and interdependencies of states across these boundaries. Complete state evolutions generated by such models may be regarded as a proxy of the real world, provided they are run at sufficiently high resolution and incorporate the most important processes. Such a simulated reality can be used to test hypotheses on the functioning of the coupled terrestrial system. Coupled simulation systems, however, face severe problems caused by the vastly different scales of the processes acting in and between the compartments of the terrestrial system, which also hinders comprehensive tests of their realism. We used the Terrestrial Systems Modeling Platform (TerrSysMP), which couples the meteorological Consortium for Small-scale Modeling (COSMO) model, the land-surface Community Land Model (CLM), and the subsurface ParFlow model, to generate a simulated catchment for a regional terrestrial system mimicking the Neckar catchment in southwest Germany, the virtual Neckar catchment. Simulations for this catchment are made for the period 2007–2015 and at a spatial resolution of 400 m for the land surface and subsurface and 1.1 km for the atmosphere. Among a discussion of modeling challenges, the model performance is evaluated based on observations covering several variables of the water cycle. We find that the simulated catchment behaves in many aspects quite close to observations of the real Neckar catchment, e.g., concerning atmospheric Published by Copernicus Publications. 4438 B. Schalge et al.: Atmosphere–land-surface–subsurface simulation dataset boundary-layer height, precipitation, and runoff. But also discrepancies become apparent, both in the ability of the model to correctly simulate some processes which still need improvement, such as overland flow, and in the realism of some observation operators like the satellite-based soil moisture sensors. The whole raw dataset is available for interested users. The dataset described here is available via the CERA database (Schalge et al., 2020): https://doi.org/10.26050/WDCC/Neckar_VCS_v1.

and Zeng et al. (2002) to better represent European tree types. Vegetation height was set to 7m 199 for needle-leaf trees and 10m for broad-leaf trees to account for partial coverage by shrubs, to 20 -120cm for crops, and to 10 200 -60 cm for grass depending on the time of the year with low values in the winter months and largest values in July and August. 201 Since we consider only one crop type, we do not specify a harvest date when the plant height drops to its minimum, but assume 202 a smooth decline between August and October. 203 For the representation of soils in CLM we use the 1:1,000,000 soil map (BUEK1000) provided by the Federal Institute for 204 Geosciences and Natural Resources -BGR 205 (http://www.bgr.bund.de/DE/Themen/Boden/Informationsgrundlagen/Bodenkundliche_Karten_Datenbanken/BUEK1000/bu 206 ek1000_node.html). This soil map is available for entire Germany; thus only small areas in Switzerland and France are missing 207 outside the Neckar catchment for which we assume a nearby soil class. BUEK1000 offers sand and clay percentages as well 208 as carbon content for two to seven soil horizons down to a maximum depth of 3m for each soil type. The carbon content is 209 used to infer soil color. For urban areas (modeled as bare soil, as mentioned above) a fixed soil color (class 8 in CLM) was 210 used. 211 Since soil properties may vary substantially at scales smaller than the 1km for which BUEK1000 is appropriate, which might (1) The BUEK1000 soil map is randomly sampled at 1995 point locations with one sample every 5 km 2 on average, a 215 minimum sample distance of 250 m, and at least one sample for each soil type of the original soil map. This strategy resulted 216 from extensive testing by minimizing the tradeoffs between reproducing the main features of the original soil map and creating 217 variability at finer resolution. 218 (2) The sample locations are used as conditional points for further interpolation. Here, texture, carbon content, and depth 219 of the first three soil horizons are extracted from the BUEK1000. In addition, the sand content of the original map was increased 220 by approximately 20% resulting in a slightly higher hydraulic conductivity because previous simulations yielded too shallow 221 unsaturated zones. 222 (3) Experimental variograms and cross-variograms are calculated for all variables and exponential models were fitted to 223 all spatial structures. 224 (4) A texture map (sand and clay percentage) is generated using a single realization based on conditional co-simulation 225 (Gomez-Hernandez and Journal, 1993) to provide the sub-scale variability (<1 km 2 ). Horizon depths and carbon content are, 226 however, assumed to have a smoothed spatial variability; therefore, they are interpolated based on ordinary kriging.
Since ParFlow describes retention and hydraulic conductivity curves based on van-Genuchten-Mualem parameters, 228 pedotransfer functions are applied to estimate these parameters. The pedotransfer functions of Cosby et al. (1984), Rawls 229 (1983) and Tóth et al. (2015) are used and selected based on data availability, applicability of the particular approaches, and 230 previous evaluations conducted in the area (Tietje and Hennings, 1996). 231 In order to keep soil porosity identical between CLM and ParFlow, we replaced the porosity calculation within CLM (which 232 uses a different pedotransfer function). The Manning's surface roughness was set to a constant value of 5.52×10 -4 h/m 1/3 and 233 the specific storage to 1×10 -3 . The chosen surface roughness value results in a realistic base flow for the local rivers without 234 calibration. Repercussions of this choice are discussed in Chapter 6. Slopes of the main rivers are additionally smoothed to 235 avoid artificial ponded areas. 236 In order to allow for realistic flow in the saturated zone, the 3-D geologic model of the geological survey of the state of Baden-237 Württemberg was used from which eleven rock types were defined for Baden-Württemberg (see Figure 4). Some characteristic 238 features of the domain, such as middle Triassic and Jurassic karst aquifers, are not included to avoid the manifold hydrological 239 challenges related to its modeling. For areas outside of Baden-Württemberg we extended the rock types at the boundary 240 outwards to cover the full computational domain. Tab. 1 summarizes porosity and hydraulic conductivity used in the domain 241 for the different stratigraphic units. As already mentioned above, karst features of limestones are not considered, and porosities 242 in stratigraphic units containing limestones and crystalline rocks are set considerably higher than in nature to somewhat counter 243 this. 244 Not covered by the discussed data sets are the large alluvial bodies filling large part of the Neckar valley throughout the domain 245 (Riva et al., 2006). Up to 30% of the runoff takes place in the subsurface especially during periods of base flow according to 246 a sub-catchment simulation performed for the year 2007. In that simulation we used measured precipitation and river discharge 247 data together with the simulated evapotranspiration to calculate the water balance over a whole year. While our simulated 248 evapotranspiration rates may be inaccurate, it is implausible that this can account for 30% of the precipitation; i.e. the water 249 could have left the domain only through the subsurface. Thus, gravel channels are needed to account for this lateral flow. Since 250 the valleys in the catchment are often small compared to the limited horizontal resolution of the model, we conceptualize the 251 alluvial bodies as gravel layers underneath all river cells (cells with a mean pressure head >0.1m) and directly next to rivers 252 (riverbanks, i.e., one grid point besides each river cell). The assumed gravel layers reach from beneath the soil down to a depth 253 of 8m. The gravel cells are parameterized with a high hydraulic conductivity of 1 m/h, a porosity of 0.6 and van-Genuchten 254 parameters of 2 for n and 4m -1 for α (residual saturation is 0.06 cm 3 /cm 3 ). Our setup results in a reasonable distribution of 255 surface and subsurface discharge at the outlet of the catchment and reasonable riveraquifer exchange fluxes. In addition to 256 the gravel channels, we included a layer of weathered bedrock, which starts below the soil and extents down to a depth of 6m. 257 This layer is characterized by substantially larger porosity (0.4) and hydraulic conductivity (0.1 m/h) than the rock below. This 258 layer was added to enhance subsurface flow and counter the common occurrence of too shallow water levels if these features 259 are not included. 260 https://doi.org/10.5194/essd-2020-24 In the following, we present example results of the virtual-reality simulations in order to demonstrate its potential for a better 265 understanding of the dynamics in coupled terrestrial systems. We will also show that the simulations quite well resemble 266 observations in the real Neckar catchment, and thus can be used to develop and evaluate modelling and prediction strategies. 267 Precipitation is the strongest hydrological driver in this region; thus its realistic spatial and temporal variability in the domain 268 including its statistical relations with topography is important. Also, the state of the atmospheric boundary layer, which reflects 269 the interaction of the land surface with the atmosphere is a critical component of the terrestrial system, which should be 270 represented by the simulation with some confidence. Along with the comparisons we will also discuss the challenges 271 experienced with such a modeling setup. 272 while July was much wetter, but the increased solar radiation and thus temperatures compared to April result in higher ET 285 rates and thus a quicker drying of the top layer of the soil. Figure 6 indicates a reduction in ET when the distance to groundwater 286 falls below 15 -100cm, depending on soil properties with faster ET reduction for increasing soil sand contents. Such relations 287 are less obvious for cells with significant plant cover: while trees show overall higher ET and almost no ET change with 288 distance to groundwater due to their deep root zones, ET variability increases with larger distances to groundwater (not shown). 289 Also crops and grassland show limited ET changes as a function of distance to groundwater, which can, however, be explained 290 by the high water availability (no water stress) in the time period considered. Figure 6 also contains a small number of grid-291 https://doi.org/10.5194/essd-2020-24 These relate most likely to cells that retain high levels of upper-level soil moisture even during dry periods to support higher 293 evaporation. This could be either because of the soil properties or due to high flow from neighboring cells. Time-series of soil 294 moisture for the 10 days preceding 30 th April and 30 th July are shown in Figure 7 to illustrate this effect. It can be seen that 295 volumetric soil moisture for cells with deep water table but high evaporation is much more similar to cells with shallow water 296 table than to cells with deep water table but low evaporation. This means that these cells can retain high upper level soil 297 moisture despite having a low groundwater level. Most of these cells are near rivers and can receive just the right amount of 298 water from lateral flow to keep evaporation high while groundwater level stays low. 299

Precipitation 300
We compare the simulated precipitation with the 1 km 2 gridded REGNIE product of DWD, derived from in-situ precipitation precipitation is dominated by advection from the west, which result in maxima over the upwind and peak zones of the 308 mountains and leeward minima. The simulated winter pattern (j) compares well with REGNIE (k), but the model 309 underestimates precipitation in the northwestern part of the catchment (l). Over the mountains a slight lateral shift of this kind 310 of precipitation pattern results in neighboring areas with under-and overestimation also found for COSMO simulations coupled 311 to its own TERRA land surface model (e.g., Dierer et al., 2009;Lindau and Simmer, 2013). In fall, the difference pattern 312 between simulations and REGNIE (i) is similar to the winter pattern, but has smaller contrasts. In spring, the simulated 313 precipitation is higher compared to REGNIE. In the summer (June to August), cloud bases are usually higher and reduce the 314 patterns caused by the luff-lee effects. Moist air extends further to the east and south and gets staunched by the alpine upland 315 leading to enhanced precipitation there. The simulated summer precipitation pattern, which is dominated by convective 316 precipitation, resembles the REGNIE pattern but exceeds the latter by about 20% lower over large parts of the catchment 317 precipitation is nearly identical between simulations and observations. The model reproduces the seasonal cycle of maximum 326 daily precipitation well, however with larger differences in the summer (see also Dierer et al. 2009). 327

Atmospheric State Variables and Surface Radiation 328
We compare the atmospheric boundary layer (ABL) of the virtual catchment to observations from the meteorological tower at 329 Karlsruhe Institute of Technology (KIT; Kalthoff and Vogel, 1992) and with DWD radiosonde observations in Stuttgart (STG) 330 (see Figure 3 for locations and Table A1 for details of observed quantities). To avoid a biased comparison due to land-cover 331 mismatches between the simulation and the actual land use at the observation sites, the simulations are averaged over five-by-332 five atmospheric grid boxes centered around the observation sites. 333 The 10m mean diurnal minimum temperatures in the virtual catchment are between 0.5 K (January) and 2.5 K (August) higher 334 than observed ( probably be further reduced by fine tuning radiation interaction parameters in CMEM, and by including orographic effects on 401 the effective incidence angle. These biases will be addressed by an improved exploitation of the uncertainty of the radiation 402 interaction parameters and by including in CMEM a two-stream approximation to better simulate cases with dense vegetation 403 in the future. 404 The microwave observations retrieved from the virtual catchment show a typical situation encountered in data assimilation; 405 more often than not there are biases between simulated and remote sensing observations. This discrepancy usually has multiple 406 causes, which can relate to the observations themselves, assumptions in the observation operator used to simulate the virtual 407 observations, and in the model used to generate the systems state variables entering the observation operator. Even if these 408 differences cannot be removed, such observations can be highly valuable for data assimilation as long as temporal tendencies 409 are meaningful information. Usually the bias is statistically corrected and thus only the information in the temporal and (if 410 meaningful) spatial variability of the observations is exploited for moving the model states towards the true states. 411

Evaluation of River Discharge 412
We compare river discharge in the virtual catchment with observations made in the Neckar catchment at the gaging stations 413 Rockenau, Lauffen, and Plochingen for a three-year period from 2007 to 2009 (Figure 16). The range of the hydrological 414 response to precipitation in the virtual catchment is in adequate agreement with the observations; this is noteworthy since no 415 calibration to runoff data has been applied to the model. The simulated discharge peaks are, however, higher and delayed by 416 about one to three days compared to the observations. A reason could be a too large Manning´s coefficient and the model 417 resolution. In the discussion we suggest a scaling of Mannings coefficient to account for the mismatch between true river width 418 and the model resolution in order to better represent realistic flood dynamics. In spring and summer, the response to 419 precipitation is significantly smoother than observed and peak amplitudes vary with respect to peak amplitudes of the 420 https://doi.org/10.5194/essd-2020-24 observations. The differences between observed and simulated precipitation discussed above and the effects of the less 421 predictable convective events during these seasons may also play a significant role. Convective events will always be displaced 422 in space and time compared to the observations and may even show different individual life cycles including lifetime and 423 amplitude. Finally, the base flow is much lower compared to the real catchment during dry periods, most likely because the 424 grid resolution is considerably larger than the actual river width and the unresolved subsurface spatial heterogeneity. An 425 increased hydraulic conductivity via an increased soil sand content may reduce the base flow further as infiltration increase s. 426 The results are further evaluated comparing the flow duration curve and the monthly run off coefficient. The former represents 427 the statistical probability to exceed a specific discharge value within a given time period while the latter is the ratio between 428 runoff and precipitation over the catchment area. Figure 17 shows the lower exceedance probability of the virtual catchment 429 compared to the observations, in particular for low discharge rates, a behavior attributed to the lower base flow component 430 and confirmed by the too low runoff coefficients in spring and summer but similar coefficients during the rest of the year 431 ( Figure 18). We hypothesize that in this period the simulation has a lower hydrological response also due to missing subsurface 432 heterogeneity. As stated above, we have neglected karst features, which are known to produce fast lateral subsurface flows. 433 Overall, the model captures the general statistical features of the catchment including the typical seasonal trends quite well, 434 while differences are noted related to hydrological extremes and base flow. These differences could be reduced by model 435 calibration from which we refrain because hydrological extremes are not primary the objective of this study. We discuss 436 options to improve the representation of river discharge further below. 437

5
Discussion 454 The size of the catchment and resolution considered (400 m) pose an enormous challenge in terms of required CPU-time. Still, 455 the applicability of Darcy's law with laboratory-based parameters can be debated as we have to resort to apparent model 456 parameters to produce realistic mass fluxes in the compartments. By compromising these technical and physical aspects in the 457 setup of the virtual catchment, we experienced several challenges; four of them will be discussed which we believe to be 458 inherent to simulating energy and mass fluxes across compartment boundaries with partial-differential-equation-based, high-459 resolution coupled models. Soil parameters: As outlined in section 2, the soil hydraulic parameters were generated based on soil maps of the real Neckar 486 catchment. According to the maps, the soils in the catchment consist mainly of clay and silt, which have rather low saturated 487 hydraulic conductivities and small air entry pressure values. In large areas of the domain, the water content in our first 488 simulations was close to saturation, even for upper soil layers, and the infiltration velocities were unrealistically low. Reasons 489 are the soil parameters, which do not capture the true soil heterogeneity; moreover, real infiltration often takes place in root 490 channels, small fractures, and other small structures. Thus, infiltration is always underpredicted by models using observed soil 491 parameters assuming homogeneity. Infiltration processes may be better captured with dual domain approaches, which are, 492 however, computationally demanding. A workaround would be to change the soil hydraulic parameters in order to obtain 493 stronger infiltration. Currently, we use an artificially increased sand percentage of the soils in order to stay consistent with the 494 concept of the pedotransfer functions used in CLM. We will also test known scaling rules (e.g., Ghanbarian et al., 2015) to 495 increase for example the saturated hydraulic conductivity for larger soil units. These rules should be applied on the soil 496 hydraulic parameters, estimated by the pedotransfer functions. 497

Conclusions and Outlook 498
In the present study we show the development and the data generated based on a integrated subsurface-land surface-atmosphere 499 system TSMP. Plausibility tests for the derived virtual reality which tries to mimic the Neckar catchment in southwestern 500 Germany, show that the virtual catchment is able to reproduce realistic behavior when compared to measurements. 501 Comparisons of simulated precipitation and ABL statistics show a very reasonable agreement with real observations. However, 502 comparisons with observed passive microwave measurements by satellites shows clearly a systematic bias which is probably 503 related to a mixture of systematic errors in the latter, assumptions in the used forward operator, parameterizations of land 504 surface properties (soil parameters) in the simulation, and missing processes therein (e.g., preferential flow, hill-slope 505 processes). The analysis also shows a realistic connection between evapotranspiration and distance to groundwater in the 506 virtual catchment, while larger deviations from reality are found for river discharge dynamics. The deficiencies could be traced 507 to the model resolution, which limits the often much smaller river widths to multiples of the model resolution, and to the way 508 river discharge is handled in the ParFlow component of TerrSysMP. study we also highlighted some limitations mainly due to the still sever technical limitation and the IT-requirements. We 518 anticipate however that more sophisticated versions of the virtual catchment (higher resolution, improved parameterization of 519 sub-scale processes as discussed above) are already in progress that could be also compared to this virtual reality in further 520 study. 521 Nevertheless, we encourage the use of this data set for investigations on data assimilation, but also the general functioning of 527 catchments including cross-environmental interactions and predictability studies can profit from such complete state evolutions 528 of the regional Earth system. 529 The TerrSysMP model is built in a modular way and users are supposed to get the component models by themselves while the 530 coupling interface is provided through a git repository (https://git.meteo.uni-bonn.de). As of now, registration is required to 531 access the TerrSysMP git and wiki page.       compared to the REGNIE data set (middle column). The difference between VR and REGNIE is shown in the right column. Figure   761