A highly modular and scale-consistent Terrestrial Systems Modeling Platform (TerrSysMP) is presented. The modeling platform consists of an atmospheric model (Consortium for Small-Scale Modeling; COSMO), a land surface model (the NCAR Community Land Model, version 3.5; CLM3.5), and a 3D variably saturated groundwater flow model (ParFlow). An external coupler (Ocean Atmosphere Sea Ice Soil, version 3.0; OASIS3) with multiple executable approaches is employed to couple the three independently developed component models, which intrinsically allows for a separation of temporal–spatial modeling scales and the coupling frequencies between the component models.
Idealized TerrSysMP simulations are presented, which focus on the interaction of key hydrologic processes, like runoff production (excess rainfall and saturation) at different hydrological modeling scales and the drawdown of the water table through groundwater pumping, with processes in the atmospheric boundary layer. The results show a strong linkage between integrated surface–groundwater dynamics, biogeophysical processes, and boundary layer evolution. The use of the mosaic approach for the hydrological component model (to resolve subgrid-scale topography) impacts simulated runoff production, soil moisture redistribution, and boundary layer evolution, which demonstrates the importance of hydrological modeling scales and thus the advantages of the coupling approach used in this study.
Real data simulations were carried out with TerrSysMP over the Rur catchment in Germany. The inclusion of the integrated surface–groundwater flow model results in systematic patterns in the root zone soil moisture, which influence exchange flux distributions and the ensuing atmospheric boundary layer development. In a first comparison to observations, the 3D model compared to the 1D model shows slightly improved predictions of surface fluxes and a strong sensitivity to the initial soil moisture content.
The evolution of the atmospheric boundary layer (ABL) is directly influenced by the spatial patterns of mass and energy fluxes from and to the land surface. Therefore, the parameterizations of land surface fluxes in atmospheric models are crucial to improving our understanding and modeling of land–atmosphere interactions, which will ultimately lead to better hydrologic predictions (e.g., Betts et al. 1996; Avissar and Pielke 1989; Chen and Avissar 1994; Sellers et al. 1997; Zeng et al. 1998; Walko et al. 2000; Chen and Dudhia 2001; Ament and Simmer 2006; Lu and Kueppers 2012). Recent studies suggest that the surface fluxes can also be strongly coupled with groundwater table dynamics (e.g., Chen and Kumar 2001; Maxwell and Miller 2005; Fan et al. 2007; Miguez-Macho et al. 2007; Niu et al. 2007; Kollet and Maxwell 2008; Maxwell and Kollet 2008; Ferguson and Maxwell 2010; Choi and Liang 2010; Niu et al. 2011). Therefore, a growing number of simulation platforms now attempt to explicitly simulate the groundwater table that is moving freely depending on the transient subsurface hydrodynamics in the coupled soil–vegetation–atmosphere system with different levels of complexity to study the linkages between land–atmosphere and subsurface hydrodynamics (e.g., Seuffert et al. 2002; York et al. 2002; Maxwell et al. 2007; Anyah et al. 2008; Jiang et al. 2009; Maxwell et al. 2011; Tian et al. 2012; Gochis et al. 2013). York et al. (2002) conducted a decadal time-scale simulation with an aquifer–soil–vegetation–atmosphere model (CLASP II-VOS-MODFLOW) to show that annually 5% (wet year) to 20% (dry year) of evapotranspiration was drawn from groundwater. For computational efficiency at decadal time scales, they used a single-column atmosphere model (~50-km grid cell) over a spatially distributed land surface hydrologic model (2-km gridcell width). Seuffert et al. (2002) coupled a regional German weather forecast model (Lokal Modell) with a land surface hydrologic model [TOPMODEL-Based Land Surface-Atmosphere Transfer Scheme (TOPLATS), based on topographic index] at 1-km resolution to demonstrate that the simulated spatial variability in soil moisture and water table depth improved predicted energy fluxes and even rain amounts. In their study, part of the improvement was also contributed by improved biogeophysics in the new land surface hydrologic model compared to the land surface scheme in the Lokal Modell. Davin et al. (2011) demonstrated improvements of simulated cloud cover, surface temperature, and precipitation at the regional climate scale, by coupling the NCAR Community Land Model version 3.5 (CLM3.5) with the Consortium for Small-Scale Modeling (COSMO) regional climate model (together referred to as COSMO-CLM). Maxwell et al. (2007) used semi-idealized simulations with a coupled groundwater–atmosphere model [Advanced Regional Prediction System (ARPS)–ParFlow] at 1-km resolution over the Little Washita catchment to demonstrate the sensitivity of thermally forced ABL development to spatially variable soil moisture and temperature. They showed that a shallower ABL is simulated over the cooler and moister river valleys, and that a certain degree of correlation exists between groundwater table depth, potential temperature in the lowest model level, and transient ABL developments. Anyah et al. (2008) used the Regional Atmospheric Modeling System (RAMS-Hydro, which included the groundwater reservoir, lateral groundwater flow, and dynamic water table river exchange) over the continental United States at seasonal time scales with grid cells of 50 km (atmosphere) and 12.5 km (land). They showed that the water-table-induced wetter soil enhanced evapotranspiration and consequently precipitation in the more arid western regions, where soil water is a limiting factor for evapotranspiration. Jiang et al. (2009) assessed the seasonal and intraseasonal evolution and pattern of precipitation over the continental United States (at 32-km resolution) with the Weather Research and Forecasting (WRF) Model coupled to a simple groundwater model (SIMGM) via the Noah land surface model (LSM). They showed that incorporating groundwater dynamics into the model produced more precipitation (1–2 mm day−1 increase in summer) over the central United States due to an increase in the latent heat flux, compared to runs without the inclusion of groundwater dynamics, which was more consistent with the observations. Maxwell et al. (2011) presented improvements in the WRF simulation platform by including a fully integrated surface–groundwater flow model. Using idealized simulations, they demonstrated that surface runoff and lateral flow change the spatial pattern of land surface fluxes. Additionally, they also used semi-idealized simulations over a small catchment to show the potential of the coupled model to develop rainfall–runoff predictions for water resource applications as well as wind-energy forecast applications. Recently, Gochis et al. (2013) presented the development of the WRF-Hydro model coupling tool, which has been designed to improve terrestrial hydrological processes in land surface models by providing a platform to couple existing groundwater models with different levels of complexity. The above studies and new developments show that the inclusion of groundwater dynamics enhances our understanding of land–atmosphere interaction through improved simulations of the spatial variability in soil moisture and groundwater table depth and, to some degree, shows the potential for improved forecasts.
The spatial variability in soil moisture and groundwater table depth also strongly depends on topography, land cover, and soil texture, which are highly heterogeneous at scales much smaller than the horizontal resolution of current mesoscale atmospheric models. The representation of this subgrid-scale heterogeneity affects grid-scale surface fluxes, along with the parameterizations of hydrological and biogeophysical processes implemented in land surface models (e.g., Schomburg et al. 2010, 2012; Bonan et al. 1993; Avissar and Pielke 1989).
In this study, we present the development and application of the modular and scale-consistent Terrestrial Systems Modeling Platform (TerrSysMP). In its current configuration, TerrSysMP comprises the numerical weather prediction model COSMO in a convection permitting configuration, COSMO-DE (Baldauf et al. 2011), CLM3.5 (Oleson et al. 2008), and the 3D variably saturated groundwater and surface water flow code ParFlow (Ashby and Falgout 1996; Jones and Woodward 2001; Kollet and Maxwell 2008). The Ocean Atmosphere Sea Ice Soil external coupler (OASIS3; Valcke 2013) is used to drive TerrSysMP and control the exchange of fluxes between each independent component model. The land-surface and groundwater model can be run in a mosaic mode at higher resolution than the atmospheric model to account for heterogeneity in land cover, soil texture, and topography.
The objective of this study is to present the major technical features of this new modeling platform and to study the effects of including integrated surface and groundwater flow in the hydrological component on land–atmosphere interactions. The goal is to demonstrate the usefulness of the scale-consistent coupling of the water and energy cycle in numerical models from the deeper subsurface, including groundwater, into the atmosphere. The paper is organized as follows. Sections 2 and 3 describe the component models and the related coupling interfaces of the modeling platform, respectively. In section 4, we present the numerical experiment designs, results, and discussions with idealized and one real-data case application including a brief discussion of the input data requirements to the system and preprocessing considerations. A summary and our conclusions are provided in section 5.
2. Model descriptions
a. COSMO atmospheric model
The COSMO model in its convection-permitting configuration, COSMO-DE (Baldauf et al. 2011), is the operational numerical weather prediction (NWP) model of the German Weather Service (Deutscher Wetterdienst; DWD), developed and maintained by an association of several European national weather services (http://www.cosmo-model.org/). The development of the model was based on an operational 7-km-resolution mesoscale weather forecast model (Lokal Modell; Steppeler et al. 2003), which is now referred to as COSMO-EU. The dynamical core of COSMO-DE uses the modified time-splitting approach of Wicker and Skamarock (2002). The equations are written in terrain-following coordinates with variable discretization using an Arakawa C grid. The physical packages (Doms et al. 2011) used in COSMO-DE consist of 1) a radiation scheme based on the one-dimensional δ-two-stream approximation of the radiative transfer equation (Ritter and Geleyn 1992); 2) a single-moment cloud microphysics scheme that predicts cloud water and ice, rainwater, snow, and graupel (Lin et al. 1983; Reinhardt and Seifert 2006); 3) a shallow-convection scheme based on the parameterization of Tiedtke (1989); 4) a level-2.5 turbulence parameterization of Mellor and Yamada (1982) and a Blackadar mixing-length scale parameterization (Blackadar 1962) resulting in a flux-gradient representation for subgrid-scale fluxes with diffusion coefficients and a turbulent length scale; 5) a surface transfer scheme within the framework of the turbulence scheme to compute the transfer coefficients for heat and momentum (Raschendorfer 2001; Doms 2001); and finally 6) a 1D multilayer soil and vegetation model TERRA (Doms et al. 2011).
b. CLM land surface model
Version 3.5 of the Community Land Model (Oleson et al. 2004, 2008) is the land surface scheme of TerrSysMP. CLM accounts for surface energy, moisture, and carbon fluxes between the shallow soil (10-layer discretization), snow (5-layer discretization), vegetation, and the atmosphere. Surface heterogeneity can be characterized with five land-cover types (i.e., the land components (glacier, lake, wetland, urban, and vegetated) within each grid cell. CLM uses a tiling approach at different subgrid levels to simulate subgrid-scale heterogeneity. Each land unit can have multiple soil columns. The third subgrid level is the plant functional type (PFT), and each soil column can have up to 4 of 16 possible PFTs (Bonan et al. 2002), each characterized by distinct plant physiological parameters. CLM3.5 prognostic variables are temperature of canopy, soil and snow layers, canopy water storage, snow depth, snow mass, snow water equivalent, and soil moisture content, among others. The biogeophysical processes simulated by CLM include absorption, reflection, and transmittance of solar radiation; absorption and emission of longwave radiation; momentum, sensible heat, and latent heat transfer from the canopy and the soil surface; plant physiology and photosynthesis; and canopy, snow, and soil hydrology.
c. ParFlow groundwater model
ParFlow is a variably saturated, groundwater flow model integrated with a two-dimensional overland flow simulator (Ashby and Falgout 1996; Kollet and Maxwell 2006). The model solves the mixed form of the 3D Richards equation and the kinematic wave equation for the overland flow boundary conditions. The Richards equation is discretized using a cell-centered finite-difference scheme and an implicit Eulerian scheme in space and time, respectively. The resulting discretized equations are solved by a Newton–Krylov nonlinear solver (Jones and Woodward 2001). The kinematic approximation of the surface flow equation is discretized using a standard upwind finite-control volume approach in space, and an implicit Eulerian scheme in time. More details on the numerical aspects and other features of the model can be found in Jones and Woodward (2001) and Kollet and Maxwell (2006).
Recently, a terrain-following grid transformation with a variable vertical discretization has been implemented in the model (Maxwell 2013). With this new feature the number of vertical layers is no longer determined by the maximum topographic height of the region to be simulated. The terrain-following grid significantly reduces the computational effort especially in cases of strong relief and allows for a finer grid resolution near the surface and within the root zone, on the order of 10−2 m, which considerably improves the simulated grid-scale fluxes at the land surface (Smirnova et al. 1997; van Dam and Feddes 2000).
3. Coupling interface
The existing coupling technologies for earth system modeling can be roughly categorized into two groups: the integrated single-executable approach and the multiple-executable approach (Valcke and Dunlap 2010). The integrated approach, which uses, for example, the Earth System Modeling Framework (ESMF; Hill et al. 2004) and the Model Coupling Toolkit (MCT; Larson et al. 2005), is more suitable in controlled development environments. The multiple-executable approach, for example, using OASIS3 (Valcke 2013), is more suitable for coupling independently developed codes, which usually experience frequent updating by their developing communities. With this in mind, the multiple-executable approach using the OASIS3 coupler was employed to develop TerrSysMP (Fig. 1), which currently consists of COSMO, CLM, and ParFlow, all of which were developed and extended independently by different research groups. TerrSysMP is highly modular and can be configured to run with different coupled component models: COSMO coupled with CLM, CLM coupled with ParFlow (using offline atmospheric forcing), and the fully coupled system (COSMO–CLM–ParFlow). Additionally, each model can be compiled and run as a stand-alone independent model within TerrSysMP. In the following section, we present the technical details of the external coupler and the coupling methodology used between the component models.
The OASIS3 coupler uses the communication techniques based on the message passing interface standards MPI1/MPI2 and the Project for Integrated Earth System Modeling (PRISM) Model Interface Library (PSMILe) for parallel communication of two-dimensional arrays between the main process (OASIS3) and the participating component models (for details, see Valcke 2013). The PSMILe library is interfaced with the existing codes of each component model; the interface consists of three main elements: 1) model initialization and definition [model grid, model partition, input/output (I/O) coupling variables], 2) sending–receiving of coupling fields, and 3) termination of simulation. For ParFlow, a new library, “oas3,” based on the existing built-in MPI1 library and the OASIS3 interface, was added. The sending and receiving of the coupling fields is done sequentially within the time stepping of the component models. The OASIS3 configuration file is used to specify the sequence of coupling, the coupling frequency, the names of the coupling fields, the spatial grid of the coupling fields, and finally the type of transformations of the 2D coupled fields.
In this study, we use the time-integration–averaging and spatial interpolation operator for the transformation of 2D variables. While the time-integration operator is only used for precipitation sent from COSMO to CLM, the time-averaging operator is used for all remaining atmospheric forcing variables. The time-integration/averaging operator is required, when the coupling frequency is greater than a component model time step. The interpolation operator is essential, when coupling variables are defined at different grid resolutions and orientations (e.g., COSMO variables are defined in a rotated geographical coordinate system and CLM variables are defined in a regular geographical coordinate system). OASIS3 uses the interpolation techniques of the Los Alamos National Laboratory Spherical Coordinate Remapping and Interpolation Package (SCRIP 1.4 library). In this study, we used the SCRIP bilinear interpolation operator for COSMO variables and the SCRIP distance-weighted averaging operator for CLM variables.
The OASIS3 interface added to COSMO and CLM exchanges the atmospheric forcing terms and the surface fluxes in a sequential manner. The atmospheric state of COSMO at its lowest level and current time step is used as the forcing term for the land surface model. The land surface model then computes the surface energy fluxes, momentum fluxes, albedo, and outgoing longwave radiation, which are sent back to COSMO (see Fig. 1). The dimensionless surface transfer coefficients of COSMO are then updated with these fluxes, and the vertical gradients at the bottom level are calculated using the surface temperature from the previous time step. The updated dimensionless surface exchange coefficients are formulated as (for a Cartesian coordinate system)
where and are dimensionless transfer coefficients for heat and momentum; SH is the sensible heat flux (W m−2), TAU is the zonal momentum flux (kg m−1 s−2); is the specific heat capacity of dry air at constant pressure; is the air density (kg m−3); and are the temperatures (K) of the surface and at the lowest model level, respectively; is the Exner function; and u and U are the zonal wind component and mean wind speed (m s−1) at the grid center, respectively. A Prandtl number of 1 is assumed, and the value of the transfer coefficient for heat is also used for the transfer coefficient for moisture.
The new surface temperature and surface humidity (qs) are estimated based on the computed outgoing longwave radiation and latent heat flux (LH), respectively, from CLM. The direct and diffuse albedo and the outgoing longwave radiation computed by CLM are also sent to COSMO as a lower boundary condition for the radiative transfer calculation. The following variables from the COSMO model are sent to CLM: air temperature (T), wind speed (U), specific humidity (QV), convective and grid-scale precipitation (Rain), pressure (P), incoming shortwave (SW) and longwave (LWdn) radiation, and measurement height. The SW radiation constitutes two components: direct and diffuse radiation.
Instead of using the existing tiling approach in CLM, we adopted the mosaic approach by defining only one PFT, one soil column, and one land unit for each CLM grid cell, which results in multiple high-resolution land surface grid cells for each atmospheric grid cell. This is important for accurately representing heterogeneity and simulating 3D subsurface hydrodynamics coupled with overland flow. This approach, however, removes the “compete for water” feature when multiple PFTs share the same soil column. We choose this explicit subgrid-scale approach because we can use the subscale pattern of the surface and subsurface for redistributing soil moisture via ParFlow (which indirectly reintroduces the compete-for-water feature).
Soil hydrology in the original CLM is 1D and only takes into account the vertical fluxes of moisture, with a simple parameterization of the groundwater table as a lower boundary condition and a decoupled surface routing along the upper boundary. However, at the catchment scale the interaction between surface and subsurface flows is a key component of the hydrological budget (e.g., Kollet and Maxwell 2006). To overcome the aforementioned limitations in the hydrological parameterization in CLM, the 3D variably saturated groundwater flow model ParFlow, with a free-surface overland flow boundary condition, is used to compute the surface runoff and subsurface hydrodynamics and thus the redistribution of soil moisture and groundwater flow in a continuum approach. The OASIS3 interface added to CLM allows dynamic interaction with the groundwater flow model and also with the atmospheric model in fully coupled mode. When coupled with ParFlow, the 1D soil column moisture prediction in CLM is replaced by the ParFlow approach (in 1D or 3D formulation). In the sequential information exchange procedure, ParFlow sends the updated relative saturation and pressure (Ψ) for the top 10 layers to CLM. In turn, CLM sends the depth-differentiated source and sink terms for soil moisture [top soil moisture flux (qrain), soil evapotranspiration (qe)] for the top 10 soil layers to ParFlow (see Fig. 1). CLM and ParFlow coupled together represent the hydrologic component of TerrSysMP.
4. Numerical experiment designs, results, and discussions
Numerical experiments with TerrSysMP for idealized and real-data cases are presented to demonstrate the system’s capability to simulate key hydrological processes, which potentially propagate to the ABL, and to examine the influence of an integrated surface–groundwater flow model on coupled mesoscale simulations. In all runs, a coupling frequency of 900 s is used between the atmospheric and hydrologic component models, which have time steps of 10 and 900 s, respectively. The coupling frequency matches the frequency of the radiation routine updates in the atmospheric model. The atmospheric variables are averaged over this time period and sent to the hydrologic component.
a. Idealized test cases
The first two idealized simulations focus on the key physical processes of runoff production (excess rainfall and saturation) at different hydrological modeling scales (with and without the mosaic approach) and illustrate how TerrSysMP redistributes soil moisture, which potentially effects ABL evolution and the local circulation. The third idealized simulation focuses on the drawdown of the water table through groundwater pumping, and its effect on land–atmosphere interactions. All simulations were performed with the fully coupled TerrSysMP system, including COSMO, CLM, and ParFlow.
1) Impact of integrated surface–groundwater flow
Two types of simulations were performed with differing hydrologic boundary conditions at the land surface. In a control simulation (CTRL RUN), integrated surface runoff and groundwater–surface water interactions were included in ParFlow via the overland flow boundary conditions. In a second simulation, surface runoff and interactions with groundwater were neglected using a simple zero flux boundary condition with the removal of ponded water at the land surface (1D RUN). Thus, in the second simulation type, the lateral distribution of moisture due to surface runoff and reinfiltration is not taken into account, which is similar to traditional land surface modeling approaches, which eventually transfer this excess water directly to a river-routing scheme.
In both numerical experiments the model geometries are identical. The domain size is 64 × 64 km2 in the horizontal with a grid spacing of 1 km and a moderate topographic slope of 3% in the negative x direction. The atmospheric domain is discretized into 64 × 64 × 50 grid cells with a vertically stretched grid and near-surface resolution of 20 m. Similarly, the subsurface domain is discretized into 64 × 64 × 30 grid cells with 10 vertically stretched layers near the surface (2–100 cm) and 20 constant levels (135 cm) extending to 30 m below the surface. In the horizontal directions, the atmosphere was initialized as homogeneous with a mean wind of U = 0.1 m s−1. The temperature was initialized to be slightly stable (T = 300 − 0.005Δz K) for z < 10 km while relative humidity was set to 50% for z < 1500 m. The land cover was specified as the CLM crop PFT (leaf area index = 1.5). The subsurface parameters were set as follows: saturated hydraulic conductivity, m h−1; van Genuchten parameters, ; n = 1.41; and porosity . The soil moisture was initialized with a hydrostatic profile perturbed in the unsaturated zone with spatially uncorrelated random noise, and a groundwater table at 2-m depth. Periodic lateral boundary conditions were used in the atmospheric model, and a no-flow lateral boundary condition was used in the subsurface model. A time step of 10 s was used for the atmospheric model. The simulation was integrated for 24 h starting at midnight. A rain rate of qr = 0.015 mm s−1 was applied in the middle of the domain over an area of 25 km2 between 0100 and 0200 LT, to generate surface runoff by infiltration excess. To ensure similar forcing for both the CTRL RUN and 1D RUN, and to remove any feedbacks associated with rainfall during the night, precipitation was prescribed by the OASIS3 coupler to the land surface model.
The CTRL RUN, which includes surface–groundwater water interactions and runoff generation due to excess infiltration, results in saturation and ponding of rainwater on the land surface, triggering overland flow. This overland flow infiltrates into the subsurface along its flow path, generating 3D soil moisture redistribution patterns in the subsurface. In the 1D RUN, ponded water is removed from the computational domain; thus, redistribution is absent.
Figure 2a shows the difference in the spatial distribution of the relative saturation ratio for the top 10 cm of soil at 1400 LT between CTRL RUN and 1D RUN. Figure 2b shows the time series of at the locations S1 (immediately downslope of the area with precipitation) and S2 (center of the area with precipitation). The soil moisture gradient in the x direction is due to the redistribution of ponded water by overland flow. The largest differences occur at S1. Similar qualitative differences in spatial patterns of soil saturation were also shown earlier by Maxwell et al. (2011), who used, however, rainfall amounts in excess of 2400 mm compared to 50 mm in this study, which is more realistic. In CTRL RUN, overland flow redistributes ponded water downslope, which cannot happen in 1D RUN (see time series of at location S1; Fig. 2b). For the grid cells, which receive rainfall (e.g., S2), the top 10 cm of soil moisture shows the same increasing trend for both simulations during the rainfall event. After the rainfall event, drops in 1D RUN, due to the removal of ponded water and the infiltration of soil moisture to deeper layers, while in CTRL RUN, remains saturated for a longer time, as the top soil layers attain a steady state between infiltration to deeper layers and surface exchange with ponded water.
Figures 2c and 2d show the vertical cross sections (y = 31 km) of the virtual potential temperature and the surface fluxes along the cross section for CTRL RUN, respectively. A cooler and moister ABL is simulated over the moister soil. The increase in soil moisture affects the partitioning of surface energy fluxes over this area, where the latent heat flux increases and the sensible heat flux decreases (Fig. 2d). Figures 2e and 2f show the differences between CTRL RUN and 1D RUN for the same cross sections. The asymmetry in surface fluxes, due to soil moisture redistribution by overland flow in CTRL RUN, is absent in 1D RUN (Fig. 2f). This asymmetry visible in difference plots for surface fluxes also propagates into the ABL (Fig. 2e). The top soil moisture for CTRL RUN is moister than for 1D RUN, which results in lower sensible heat fluxes (SHFs) and higher latent heat fluxes (LHFs) compared to those in 1D RUN. Especially at S1, where overland flow redistributes the ponded water from uphill grid boxes, SHF is lower. Accordingly, the simulated ABL at 1400 LT is slightly cooler for CTRL RUN compared to 1D RUN, with the presence of a slightly warmer entrainment zone above the ABL (see Fig. 2e).
In both runs, the juxtaposition of wet and dry patches produces a local circulation, which intensifies at 1400 LST. The results are similar to circulations created by introducing a patch of unstressed vegetation into a region with bare dry soil, as shown in the idealized 2D simulations of Avissar and Pielke (1989). Here, we want to emphasize the differences in the magnitude of the circulation and the driving fluxes between the two simulations. Stronger downdrafts are simulated for CTRL RUN (Fig. 2g) due to lower sensible heating from the moister soil as a result of overland flow, reinfiltration, and soil moisture redistribution. The asymmetry in the vertical cross section of soil moisture is also clearly visible in the vertical cross section of the wind vectors. Thus, the spatial gradient of the soil moisture simulated by an improved representation of infiltration and the 3D soil moisture redistribution also affects the circulation and ABL evolution.
2) Effect of subgrid-scale topography
The goal of this test case is to demonstrate the importance of subgrid-scale, topographically driven surface-subsurface flow on land–atmosphere interactions. The local topographic slopes are one of the key parameters that dictate the time scales of overland flow and groundwater flow. Commonly, due to computational costs, coarser resolutions are adapted depending on the problem size. This coarsening of model scales (grid-cell size) also leads to smoothing of the local slopes, resulting in moister domains (e.g., Kuo et al. 1999; Zhang and Montgomery 1994; Sulis et al. 2011). Here, we investigate how this smoothing affects the soil moisture distribution and land–atmosphere interactions. We use CTRL RUN presented in section 4a(1) as a reference case, and compare with our SUBGRID RUN. We use a mosaic approach in SUBGRID RUN by changing the model scale (grid size) of the hydrologic component. The model scale is disaggregated from 1 km to 200 m (domain size of 320 × 320) with a heterogeneous slope pattern mimicking subgrid-scale channels or rivulets (Fig. 3a). The mean slope of these high-resolution grid cells when aggregated to 1-km resolution conserves the magnitude and direction of slopes as in CTRL RUN. The results of SUBGRID RUN for the hydrologic component were aggregated to 1-km resolution for comparison with CTRL RUN. In SUBGRID RUN, overland flow extends farther outward from the area with rainfall (up to x = 24 km) and redistributes soil moisture in the upper layers, while for CTRL RUN, overland flow barely exceeds the adjacent grid cells (x = 28 km). This is clearly illustrated by the difference in the top 10 cm of soil moisture between SUBGRID RUN and CTRL RUN (Fig. 3b). The time series of the top 10 cm of soil moisture in the center (S2), immediately downslope (S1), and 3 km farther downslope of the rainy area (S3) also demonstrates the reduction in the time scale of ponding and the faster redistribution of ponded water caused by the convergence pattern following the topographic slopes (Fig. 3c). For example, at S2, the peak in saturation occurs 4 h earlier in SUBGRID RUN compared to CTRL RUN. Similarly, the amplitude in soil moisture at S1 is much lower for SUBGRID RUN. Thus, the slope heterogeneity enhances the overland flow, thereby reducing infiltration and mean soil moisture content. The simulated fluxes and the ABL state react to the inclusion of the subgrid (relative to the atmospheric model) surface and subsurface heterogeneity (Figs. 3d,e). The downslope extension of the moister area in SUBGRID RUN lowers the sensible heat flux also in this area (24 km < x > 28 km), resulting in a cooler ABL and a warmer entrainment zone, as well as in a wider extent of the downdraft area (see difference plots of wind vectors and vertical soil moisture profiles between SUBGRID RUN and CTRL RUN at 1400 LT in Fig. 3f). The results from this idealized numerical experiment illustrate that subgrid-scale topography can have significant influence on the redistribution of soil moisture, which directly affects the surface energy fluxes, the local circulation, and ABL evolution, and hence supports the usefulness of fully coupled mosaic simulations.
3) Dynamic groundwater table
The impact of subsurface hydrodynamics on the lower atmosphere at daily time scales is demonstrated by an idealized simulation of groundwater abstraction, which results in a local drawdown of the free groundwater table at the kilometer scale. The size of the flat domain is 64 × 64 km2 in the horizontal with grid spacing of 1 km. The subsurface is vertically resolved in 10 stretched layers near the surface (2–100 cm) and 5 constant-depth levels (240 cm) extending down to 15 m below the land surface. For the atmosphere, discretization and initialization configurations from the previous test case were used. A uniform CLM crop PFT (leaf area index = 1.5) over a homogeneous sandy soil was assumed for the land surface and subsurface model with porosity , van Genuchten parameters , and saturated hydraulic conductivity . Soil moisture was initialized with a hydrostatic profile perturbed in the unsaturated zone with a spatially uncorrelated random noise (, for pressure head) and a groundwater table at 1-m depth. Water-level drawdown was generated by the simulation of four pumping wells, screened at the bottom of the saturated zone (aquifer), which pump at a total rate of 222 L s−1 over an area of 25 km2. This test case reflects conditions commonly observed in heavily irrigated agricultural settings like those in the U.S. Great Plains, where groundwater pumping for irrigation has resulted in water-table drops of 3–60 m below pre-irrigation levels (Sophocleous 2000; York et al. 2002). In practice, however, the pumped water is often applied directly at the land surface for irrigation or used for other purposes, while in this idealized simulation it is removed from the domain.
At 1 km × 1 km spatial grid resolution, the groundwater pumping results in a drawdown of the groundwater table of approximately 0.65 m, which leads to a fast drying of the soil in the overlying unsaturated zone (Fig. 4a). As expected, a steeper drop in soil moisture is observed in locations within the area of groundwater pumping (S4) compared to locations outside the region (S3), where reactions to this disturbance are much slower (Fig. 4b). The result is a relatively dry patch of soil, which causes a differential partitioning of land surface fluxes. Figures 4c and 4d show the anomalies of virtual potential temperature and surface fluxes along the cross section at y = 31 km at 1400 LT. The anomalies are calculated based on the mean values of the virtual potential temperature and fluxes for the given cross section at 1400 LT. A larger (by almost 20 W m−2) sensible heat flux is simulated over the drier patch due to the faster heating rate of the dry soil, which produces a warmer ABL (+0.08 K on average) above the patch. The dry patch extends down to 1 m below the land surface, because of the lowering water table (Fig. 4e) due to pumping. The differential heating of the lower ABL produces temperature gradients, which generate a local circulation pattern. Near the surface, the flow is directed from the wet to the dry soil and creates a convergence at the center of the domain. This circulation pattern is already well developed at 1400 LT (Fig. 4e) and suggests that groundwater-table drawdown (e.g., caused by groundwater pumping) may affect local circulation patterns and ABL evolution.
b. Real-data test case
The goal of this numerical experiment is to examine the effect of coupling the water and energy cycle in numerical models from the deeper subsurface, including groundwater into the atmosphere at mesoscale catchment scales, and to show the usefulness of the new simulation platform TerrSysMP in improving our understanding of land–atmosphere interactions under realistic conditions.
1) Experiment design
The numerical simulation was carried out for a week from 21 to 27 July 2012. Clear-sky conditions without rainfall over most of the model domain during that period are well suited to investigating the influence of topography-induced subsurface hydrodynamics on soil moisture distribution, fluxes, and the ABL. The model domain encompasses the Rur catchment [study region of TR32; see, e.g., Vereecken et al. (2010)], which is located mainly in Germany along its borders with Belgium, the Netherlands, and Luxembourg (Fig. 5). Surface fluxes are examined at four different stations within the catchment: 1) Selhausen (50.8658°N, 6.45°E), 2) Rollesbroich site I (50.62°N, 6.30°E), 3) Rollesbroich site II (50.62°N, 6.30°E), and 4) Merzenhausen (50.93°N, 6.29°E) and also compared with measured fluxes at these locations.
The numerical domain covers an area of 150 km × 150 km with the atmospheric domain horizontally resolved in 1 km × 1 km grid columns. Fifty levels were used in the vertical with a near-surface-layer depth of 20 m. The hydrological domain was discretized with a higher horizontal resolution of 500 m in order to better represent the heterogeneity based on the land-use data from the Moderate Resolution Imaging Spectroradiometer (MODIS) and the topography. This resolution is still coarse for the hydrological component model, and we acknowledge this limitation in our current study. We argue, however, that this scale separation is already sufficient to demonstrate the usefulness of the mosaic approach for physically based surface–groundwater flow models as shown earlier using the idealized numerical experiments. In the vertical, 30 levels were used for the hydrologic component with 10 stretched layers near the surface (2–100 cm) and 20 constant levels (135 cm) extending to 30 m below the surface.
The MODIS land cover type 5 (PFT classification scheme) provided by the Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC) was used to create the PFT map for the model domain (Fig. 6a). Urban areas were replaced with crop land cover for this study, with a leaf area index of 0.62, which is consistent with the land cover used in the operational COSMO version. The hydrological delineation of the Rur catchment is also shown in Fig. 6a. Small patches of forested land with needleleaf trees (nle) and broadleaf deciduous trees (bld) are embedded in larger areas of crops (c1) and deciduous broadleaf shrubs (c3g). The hilly landscape of the Eifel low-mountain range dominates the southern part of the model domain while the northern part is characterized by predominantly flat terrain. Mountainous areas are mostly covered with forest whereas the valleys are dominated by crops. In addition to the PFT classification, CLM also requires predefined distributions of leaf area index (LAI) and stem area index (SAI). A phenology study was conducted for the dominant vegetation types from 2002 to 2011 using the cloud-screened MODIS (MCD15A2) 8-day-composite 1-km LAI from the Aqua and Terra satellites. The spatially and temporally averaged monthly LAIs of the PFTs were used to specify the monthly leaf area index for each PFT (Figs. 6b,c). The SAI was estimated based on LAI following Zeng et al. (2002) and Sellers et al. (1996). For the sand and clay percentages, the soil texture map of COSMO was used, which is based on the Digital Soil Map of the World [Food and Agriculture Organization of the United Nations (FAO); Fig. 6d]. The soil parameters used in this study are described in Table 1. This soil texture was used for the upper 10 soil layers, while for the layers below, the soil texture information from Gleeson et al. (2011) was used. In addition to soil texture information, ParFlow needs the local topographic slopes in x and y directions. The r.watershed tool in the Geographic Resources Analysis Support System (GRASS; Neteler et al. 2012) was combined with the NCAR Command Language (NCL; NCAR 2013) to generate a preprocessing tool for creating the slopes for the given digital elevation model.
Initial and lateral boundary conditions for the atmospheric model were obtained from the COSMO-DE analysis files provided by DWD. Initial soil moisture and soil temperature were obtained from a spinup using the hydrologic component of TerrSysMP. Three sets of simulations were conducted: 1) using soil moisture–temperature from a 3D spinup (BASE RUN), 2) using soil moisture–temperature from a 1D spinup with the groundwater-table depth at 3 m (GWT-3m), and 3) using soil moisture–temperature from a 1D spinup with the groundwater-table depth at 1 m (GWT-1m). The 1D spinup was carried out by using the hydrologic component as a 1D column model by setting zero slopes, no lateral flow, the prescribed groundwater-table depth as a lower boundary condition, no overland flow, and the removal of ponded water at the land surface.
The spatial distribution of the topographic index for the Rur catchment at 500-m resolution (Fig. 7a, which also includes the location of the stations used for this study) illustrates the river network and the gradient from the hilly landscape of the Eifel low-mountain range in the south to the more flat northern part of the catchment. The simulated volumetric soil moisture of the top 10 cm of the soil layer for BASE RUN, GWT-3m, and GWT-1m are shown in Figs. 7b–d, respectively (after 40 h of model integration time). While the spatial pattern of soil moisture for BASE RUN is clearly correlated with the topographic index, this topographic control is completely absent for the two runs with a fixed water table and no lateral surface and subsurface flow. In these runs, the gridcell surface soil moisture is to some degree controlled by the fixed location of the groundwater table. At the catchment scale, the spatial mean soil moisture values for the three runs are 0.386, 0.258, and 0.338 m3 m−3, respectively (Figs. 7b–d). As expected, GWT-3m results in drier soil moisture by almost 33% compared to the BASE RUN. The GWT-1m run simulates only slightly drier soil moisture than BASE RUN (~5%) at the catchment scale. This shows the relative importance of the lower boundary condition for the simulation of soil moisture in 1D column models, which is explicitly simulated in the BASE RUN. It also clearly shows that the inclusion of a 3D groundwater model contributes more water to the dry soil through 3D redistribution, which is consistent with results from earlier studies (e.g., Seuffert et al. 2002; Maxwell et al. 2007; Jiang et al. 2009).
The simulated surface fluxes for the three runs are compared to observations at four different stations (Fig. 8). At Selhausen (Fig. 8a), BASE RUN better captures the diurnal cycle of the observed latent heat and sensible heat fluxes compared to the other two runs. At this location, the surface soil moisture results (top 10 cm) for GWT-3m and GWT-1m are drier than for BASE RUN by 30% and 20%, respectively. Both GWT-3m and GWT-1m simulate higher sensible heat fluxes and much lower latent heat fluxes. Previously, Chen and Dudhia (2001) also showed that a 10% change in soil moisture results in a 30 W m−2 variation in surface fluxes. In their simulation, they suggested that the influence of initial soil moisture was carried over into the 24–48-h simulation period. Here, the effect appears to be persistent over 6 days. At the two Rollesbroich sites (Figs. 8b,c), BASE RUN and the GWT-1m run simulate similar fluxes due to small differences in their initial soil moisture settings (5%). Both runs simulate the observed latent heat flux well; however, some of the daytime peaks are underestimated. Sensible heat fluxes are, however, overestimated, which is probably caused by an inappropriate choice of vegetation type. Currently, the land surface model has only prescribed crop parameters for wheat while the agricultural crops in the region can be barley, sugar beet, potatoes, and others. The initial soil moisture for GWT-3m is again 29% drier than for BASE RUN and underestimates the latent heat flux and overestimates the sensible heat flux. At Merzenhausen (Fig. 8d), GWT-3m simulates the observed fluxes better than do the other two simulations. The initial soil moisture results for GWT-3m and GWT-1m are, respectively, 30% and 20% drier than for the BASE RUN. BASE RUN overestimates the latent heat flux and underestimates the sensible heat flux due to higher initial soil moisture, which is potentially due to the slower-than-real drainage. A drainage reduction can be expected, when the spatial resolution of the model does not capture the true variability of the slope distribution, as discussed in the idealized test case (SUBGRID RUN). In general for the four stations, BASE RUN simulates higher latent heat flux compared to the other two runs due to the increase in root-zone soil moisture with the inclusion of the 3D groundwater model.
Finally, the variability in the diurnal cycle of the simulated ABL between the three runs at the four different stations is also examined. Figures 9a and 9b compare the diurnal cycles of the simulated potential temperature and vapor mixing ratio (6-day average) between BASE RUN and GWT-3m and GWT-1m, respectively. BASE RUN simulates a cooler and moister boundary layer compared to the other two runs. Differences are obvious from the initial growth of the ABL at sunrise until its collapse in the evening. BASE RUN simulates a relatively cooler and moister ABL, with a warm and dry entrainment zone on top. The ABL potential temperature is on average 0.5 (0.2) K lower in BASE RUN compared to GWT-3m (GWT-1m) due to the higher soil moisture and ensuing lower sensible heat flux simulated by BASE RUN, which is consistent. Similarly, the ABL vapor mixing ratio is on average 0.35 (0.2) g kg−1 higher in BASE RUN compared to GWT-3m (GWT-1m) due to higher latent heat flux. At Rollesbroich, where the initial soil moisture was different only by 5%, the difference in the amplitude of the virtual potential temperature is much weaker compared to other stations where the difference in the initial soil moisture is on the order of 20% (Fig. 9b).
The real-data case study suggests that the inclusion of a 3D groundwater model enhances the redistribution of moisture to the dry soil, which increases the overall latent heat flux, and thereby produces a cooler and moister ABL. The findings of this study are limited to short-time-scale coupled simulations, because our emphasis was on presenting the feasibility and usefulness of this system in studying the impacts of including a 3D groundwater model on land–atmosphere interactions. A longer-term (seasonal scale) coupled simulation is needed to further explore and deepen our understanding of the effect of the spatial–temporal variability of groundwater-table depth on land–atmosphere interactions.
5. Summary and conclusions
This paper presents the formulation and setup of the Terrestrial Systems Modeling Platform (TerrSysMP), which consists of the atmospheric model COSMO, the land surface model CLM, and the 3D variably saturated groundwater flow model ParFlow, coupled using the multiple-executable approach with the external coupler OASIS3. TerrSysMP is highly modular by design, and can be a very valuable tool for studying land–atmosphere interactions with explicit linkages to groundwater dynamics. The use of the OASIS3 coupler also allows for an explicit separation of spatial–temporal modeling scales and coupling frequency between the component models, allowing, for example, the land surface and groundwater model to be run at much higher spatial resolution (mosaic approach) and larger time steps compared to the atmospheric model.
Idealized simulations presented in this study focus on the key physical processes of runoff production (excess rainfall and saturation) at different hydrological modeling scales and the drawdown of the water table through groundwater pumping. The results from the coupled simulations demonstrate the propagation of these physical processes to the atmospheric boundary layer, indicating a strong linkage between integrated surface–groundwater dynamics, biogeophysical processes, and ABL evolution. The use of a mosaic approach for the hydrological component model (to resolve the subgrid-scale topography) strongly influences the simulated runoff production. While the distribution of the subgrid-scale slope may vary spatially, we argue that the subgrid-scale slopes exert a certain degree of control on the simulated mean grid soil moisture. In this study, the use of the mosaic approach affected the mean grid-cell soil moisture redistribution to a wider extent and consequently also resulted in a wider extent of a cooler ABL and a warmer entrainment zone. These results show the importance of hydrological modeling scales and also the advantages of the mosaic approach used in this study.
The impact of including a dynamic groundwater model on the surface energy fluxes and ABL evolution on the catchment scale is examined based on a real-data case simulation. The 3D hydrologic component model produces a topographically influenced soil moisture distribution, with a higher average soil moisture content compared to the 1D hydrologic component model at catchment scale. This difference in average soil moisture content is due to the availability of groundwater to draw upon, illustrating that the dynamic groundwater-table depth exerts a certain degree of control over the root zone soil moisture. The simulated surface fluxes were also compared with measured fluxes at the land surface. In general, the fully coupled reference run showed an improved level of prediction by the simulated surface fluxes, except at one location where the latent heat flux was overestimated and the sensible heat flux was underestimated. This could probably be due to higher initial soil moisture resulting from drainage reduction when the spatial resolution of the model does not capture the true variability of the slope distribution. This limitation has been acknowledged for the real-data case simulation presented in this study. The results were also found to be sensitive to the initial soil moisture content obtained from the spinup. When the root zone soil moisture difference was within 5%, the coupled runs with a 3D and 1D hydrological component model produced similar diurnal cycles of surface fluxes for the short-term simulation. This also suggests the importance of groundwater-table depth as a lower boundary condition, which is explicitly simulated by the 3D hydrological component model. Kollet and Maxwell (2008) showed that the surface fluxes are sensitive to groundwater dynamics, when the groundwater level is in the range of a critical depth that extends from 1 to 5 m. The change in the partitioning of surface energy fluxes also contributes to a relatively cooler and moister boundary layer for the coupled run with the 3D hydrologic component model, showing how the groundwater table exerts control over the ABL’s evolution. Since the circulation patterns in the ABL, which may trigger convection and precipitation, are initiated and/or moderated by the surface fluxes, a more realistic representation of the spatial–temporal variability of the soil moisture distribution is important for a better understanding and eventually for the improved prediction of the state of the coupled terrestrial system. Our results suggest that TerrSysMP is a useful tool for such studies.
Future work includes the investigation and validation of TerrSysMP with seasonal-scale simulations using a high-resolution hydrological component (<100 m) and to examine the impacts of incorporating groundwater dynamics on simulated precipitation. Currently, tests are carried out with the addition of new plant functional types in the land surface model to address land-cover heterogeneity, especially with regard to agricultural crops. Also, the land–atmosphere CO2 exchange is being implemented into the modeling platform and tested. Additionally, the coupling routines have been upgraded to use OASIS3-MCT (Valcke et al. 2013) for faster data exchange between the component models, and tests are being carried out for weak-scaling studies.
The development of the TerrSysMP was performed within SFB/TR32 (www.tr32.de), “Patterns in Soil–Vegetation–Atmosphere Systems: Monitoring, Modeling, and Data-Assimilation,” funded by the Deutsche Forschungsgemeinschaft (DFG). We thank the Deutscher Wetterdienst for the COSMO analysis data and the COSMO model code. The CLM source codes were obtained from NCAR (http://www.cgd.ucar.edu/tss/clm/distribution/clm3.5/index.html). We are also thankful to Reed Maxwell for providing us with the terrain-following grid version of ParFlow with variable vertical discretization. Finally, we would also like to thank Eric Maissonave and Edouard Davin for their initial support on the OASIS interface. The data analysis, including the pre- and postprocessing of input data for TerrSysMP, was done using the NCAR Command language (version 6.1.2).