The boundary layer, land surface, and subsurface are important coevolving components of hydrologic systems. While previous studies have examined the connections between soil moisture, groundwater, and the atmosphere, the atmospheric response to regional water-table drawdown has received less attention. To address this question, a coupled hydrologic–atmospheric model [Parallel Flow hydrologic model (ParFlow) and WRF] was used to simulate the San Joaquin River watershed of central California. This study focuses specifically on the planetary boundary layer (PBL) in simulations with two imposed water-table configurations: a high water table mimicking natural conditions and a lowered water table reflecting historic groundwater extraction in California’s Central Valley, although effect of irrigation was not simulated. An ensemble of simulations including three boundary layer schemes and six initial conditions was performed for both water-table conditions to assess conceptual and initial condition uncertainty. Results show that increased regional water-table depth is associated with a significant increase in peak PBL height for both initial condition and boundary layer scheme conditions, although the choice of scheme interacts to affect the magnitude of peak PBL height change. Analysis of simulated land surface fluxes shows the change in PBL height can be attributed to decreasing midday evaporative fraction under lowered water-table conditions. Furthermore, the sensitivity of PBL height to changes in water-table depth appears to depend on local water-table variation within 10 m of the land surface and the regional average water-table depth. Finally, soil moisture changes associated with lowered water tables are linked to changes in PBL circulation as indicated by vertical winds and turbulence kinetic energy.
The atmosphere and terrestrial hydrosphere continually coevolve through their interactions at the land surface. As the lower boundary of the atmosphere, the planetary boundary layer (PBL) affects heat, moisture, and momentum fluxes and is an important factor controlling land–atmosphere interactions and much of the weather we experience (Stull 1988). Studies of the PBL in the context of land–atmosphere interaction are often framed as involving near-surface soil moisture and atmospheric responses like temperature, boundary layer development, and precipitation. While a growing body of literature has addressed many of the complexities associated with such land–atmosphere coupling through soil moisture (Findell and Eltahir 1997; Patton et al. 2005; Juang et al. 2007; Adler et al. 2011), the role of the deeper subsurface has received less attention. Similarly, while the role of groundwater flow systems in maintaining spatial and temporal structure of soil moisture and associated surface energy partitioning has been well explored (Chen and Hu 2004; Miguez-Macho et al. 2007; Kollet and Maxwell 2008; Maxwell and Kollet 2008; Rihani et al. 2010), significant questions remain regarding how such interactions may be influenced by dynamic atmospheric feedbacks.
A portion of the previous research examining deeper subsurface–atmosphere interactions has been done using simplified representations of the hydrology (Seuffert et al. 2002; Jiang et al. 2009), the atmosphere (Bonetti et al. 2015; York et al. 2002), or both (Quinn et al. 1995). A general finding from such studies has been that extending the treatment of groundwater in coupled atmospheric models can improve terrestrial hydrologic dynamics in the form of a more capacious store and buffer of atmospheric water (sustaining streamflow and evapotranspiration) and as a means of distributing and maintaining soil moisture more appropriately across the modeled landscape. Although constrained in some aspects by the limited physics used, the insights gained from such simplified simulations serve as useful comparisons to more process-intense modeling studies. Of specific relevance for this study, Quinn et al. (1995) and Bonetti et al. (2015) demonstrated a direct correlation between water-table depth (WTD) and boundary layer height over a range of effective water-table depths on the order of several meters, despite using a simplified slab boundary layer model.
Sophisticated hydrologic and atmospheric models are increasingly being used together to further explore the pathways that link groundwater, soil moisture, and atmospheric behavior. Such models have been used to initialize land surface–atmospheric models with soil moisture profiles consistent with the regional groundwater and surface hydrology (i.e., historical climate, terrain, vegetation, etc.). For example, Rihani et al. (2015) used this approach for an idealized domain simulated with the Advanced Regional Prediction System (ARPS) to reveal the transient importance of topography, soil moisture, and water-table depth on boundary layer development, noting a correlation between PBL height and water-table depth extending to 10 m. Rahman et al. (2015) and Shrestha et al. (2014) compared atmospheric response in coupled Parallel Flow hydrologic model (ParFlow)–CLM–COSMO simulations of real-world domains that alternately reflected three-dimensional subsurface flow and limited one-dimensional flow. They demonstrated the importance of lateral surface and subsurface flow for mediation of latent heat flux in time and space and improving matches to surface energy flux measurements. Similarly, Maxwell et al. (2007) used the coupled ParFlow–ARPS model to show that surface fluxes, states, and boundary layer development are sensitive to hydrologically consistent soil moisture initializations, but that the effect of such initial consistency may decay over extended simulation (~36 h) if not maintained by lateral flow mechanisms. Finally, Keune et al. (2016) show, through a series of Terrestrial Systems Modeling Platform (TerrSysMP) simulations over continental Europe, that different subsurface parameter and process representations affect land surface fluxes as well as the coupling and feedback processes between the subsurface and atmosphere, noting a complex interaction of sensitivities to the presence of lateral flow and hydraulic parameter values.
The growing evidence for connection between groundwater and the atmosphere suggests that anthropogenic impacts to groundwater systems may extend beyond terrestrial hydrology. Extensive groundwater extraction to support agriculture, industry, and domestic use has produced regions of considerable aquifer depletion around the world (Döll et al. 2012; Scanlon et al. 2012; Wada et al. 2010). In particular, groundwater extraction to support irrigation has been linked with significant declines in important regional aquifer systems (Wada et al. 2012), with California’s Central Valley being a notable example in the United States (Famiglietti et al. 2011; Scanlon et al. 2012). Although anomalous soil moisture associated with irrigation in arid and semiarid regions has been shown to propagate effects to the atmosphere (Lobell et al. 2009; DeAngelis et al. 2010; Kueppers and Snyder 2012; Lo and Famiglietti 2013), the component effect of water-table drawdown for the same areas has received little examination and provides a key motivation for this study. Understanding the role of the water table in such regions has particular relevance, as irrigation may be periodically limited as a result of drought, water-table drawdown beyond the reach of wells, changes in regulation, or some combination of factors.
We use the subsurface–land surface–atmospheric modeling framework that couples ParFlow and the Weather Research and Forecasting Model (PF.WRF; Maxwell et al. 2011; Williams and Maxwell 2011) to analyze the impact of a lowered water table on the planetary boundary layer over the San Joaquin River basin in central California. The approach used in this study comprises a series of simulation components, described in the following sections.
a. Domain and time period of interest
The model domain is centered on the San Joaquin River basin and covers the majority of the contributing area of the watershed in the Sierra Nevada and the Central Valley floor. The domain is 270 km in the east–west dimension and 240 km in the north–south dimension. The map in Fig. 1 shows the location and configuration of the domain. The lateral discretization of the model domain for the coupled simulations is a uniform 1 km in the x and y dimensions. Vertical discretization differs for the ParFlow and WRF components, as described below. Modern vegetative cover in the region is dominated by irrigated cropland, pasture, savannah, and shrubland on the valley floor and edges that grades into evergreen needleleaf and mixed forests in the middle elevations of the Sierra Nevada and sparse or barren cover at higher elevations (NASA 2001). This land-cover distribution is reflected in the vegetation types assigned to the Noah land surface model used in the coupled PF.WRF simulations.
We selected a 2-week period in July 2003 (14–28 July 2003) as the period of interest for the coupled PF.WRF simulations. We chose to simulate a 2-week span because coupled simulation times were short enough to make completing an ensemble of runs feasible while allowing the simulation to capture a wider range of day-to-day variability. The specific 2-week span was chosen for several reasons. First, this midsummer time period is sufficiently distant from peak snowmelt season to reduce the intensity of surface flows discharging from the Sierra Nevada, flows that might otherwise dominate the Central Valley hydrologic regime. Second, midsummer in a Mediterranean climate such as California’s tends to be characterized by fairly consistent atmospheric conditions with precipitation in the mountains or valley driven by local convection rather than dominant synoptic-scale events, for example, influx of large, moist Pacific air masses. Third, the 14 days spanning 14–28 July 2003, while mostly dry, do include some small precipitation events across the Sierra Nevada crest and even extending into portions of the Central Valley, as indicated by station and radar records. These events introduce some variability into the expected atmospheric structure while the 2-week simulation period provides an opportunity for coevolution of the subsurface–land surface–atmospheric system.
b. WRF nests
The 1-km lateral-resolution PF.WRF domain is forced with lateral boundary conditions derived through offline downscaling of WRF, version 3.3 (Skamarock et al. 2008), simulations using three nests at increasing resolution (27, 9, and 3 km). The downscaling runs were done using WRF and the Noah land surface model but without the ParFlow simulation component. The extents and locations of these nests are shown in Fig. 2, with the 1-km PF.WRF domain shown with a thicker outline in the center. The WRF domains were discretized into 65 vertical layers with hyperbolic tangential stretching applied to concentrate more layers in the lower portion of the atmospheric column. Using this scheme, the lowest 2 km was discretized into 25 layers. The North American Regional Reanalysis (NARR) product (Mesinger et al. 2006; NCEP 2005) was used to drive the outermost (27 km resolution) nest over the period from 1 June through 1 August 2003. One-way nesting was employed for the interior domains. The innermost (3 km) nest was run in convection-permitting mode while convection parameterizations were employed for the coarser 27- and 9-km nests. The choice of physics schemes and other input options are summarized in Table 1. The numbers in brackets correspond to the options used in the WRF input file.
For each nested WRF simulation, a spatial and temporal (14–28 July 2003) subset of the simulated conditions from the 3-km nest was processed to produce boundary and initial conditions for the coupled San Joaquin basin domain.
c. ParFlow hydrologic model
The ParFlow component of this study simulates variably saturated subsurface flow and overland flow with a fully implicit solution of the three-dimensional mixed form of the Richards equation and a kinematic wave approximation for surface flow (Ashby and Falgout 1996; Jones and Woodward 2001; Kollet and Maxwell 2006; Maxwell 2013). The terrestrial domain is discretized at a 1-km lateral resolution and a variable resolution over a 500-m vertical depth using five layers. The top four ParFlow layers are discretized according to the scheme used in the Noah land surface model (Ek et al. 2003, p. 200): 0.1, 0.3, 0.6, and 1.0 m, in order of increasing depth from the land surface. The fifth (bottom) layer spans 498 m and functions as a vertically integrated approximation of deeper subsurface flow. A terrain-following transform (Maxwell 2013) is applied to vertical dimension of the grid such that the topographic variability of the landscape is represented across all five model layers. This conceptualization of the subsurface flow system does not represent variation in vertical gradients in depth below 2 m but captures the large-scale horizontal gradients responsible for lateral transport of water across a regional subsurface flow system like the San Joaquin River basin.
The San Joaquin ParFlow domain is defined with no-flow conditions on the lateral and bottom boundaries. An overland flow boundary condition is assigned to the top boundary, coincident with the modeled land surface. Overland flow routing is simulated by the kinematic wave approximation (Kollet and Maxwell 2006) with inputs comprising surface slopes in the x and y directions and Manning’s roughness values. Following the method of Barnes et al. (2016), slopes are derived from a ⅓-arc-s (~10 m) DEM aggregated to the 1-km-resolution grid using the watershed analysis tool in the Geographic Resources Analysis Support System (GRASS) GIS (Gesch 2007; Gesch et al. 2002; Ehlschlaeger and Metz 2014) and Manning’s roughness values assigned from reference values (Chow 2009) based on general vegetation type and topographic region.
Subsurface properties are assigned according to indicator categories that describe a hydrostratigraphic model of the system. The simplified hydrostratigraphy used here is based on datasets developed for previous studies of the region (Faunt 2009; Mansoor 2009) and includes detailed surface layers and an aggregated deeper subsurface layer that is the mode of hydrostratigraphic units over that interval. Hydraulic parameter values (permeability, van Genuchten constitutive relationship parameters, porosity, and specific storage) are assigned according to the previous studies where available, that is, near-surface layers and for the Central Valley, and using reference values for the portion of the subsurface representing the Sierra Nevada and Coastal Range mountain blocks.
The ParFlow domain is initialized for use in coupled PF.WRF simulations through a spinup process whereby surface and subsurface hydrology is dynamically equilibrated with a given climatological forcing. We initially applied a constant (in time) precipitation field to the surface of the dry domain and allowed the subsurface to wet until a regional water table and surface features began to form. We continued the spinup process by iteratively simulating the same year using an average, but transient hourly, atmospheric forcing with the CLM-based land surface model component activated (Kollet and Maxwell 2008; Maxwell and Miller 2005) until the terrestrial hydrology approached dynamic equilibrium with the forcing, as measured by the annual change in total surface and subsurface storage following documented metrics for integrated model spinup (Ajami et al. 2014a,b; Seck et al. 2015). For this study, we applied the spinup process until the change in total storage between the end of a spinup year and the beginning was less than 0.1% of the beginning storage. This criterion was applied to both hydrology scenarios described in the next section.
Terrestrial hydrology scenarios for PF.WRF
We developed two terrestrial hydrologic initial conditions for use with the PF.WRF simulations: a predevelopment (natural or no pump) condition and a modern water-table drawdown (pumped) condition. The first represents a system driven by modern meteorology but lacking local anthropogenic hydrologic influences like surface water impoundments, channel diversions, canals, groundwater extraction, and irrigation. The second represents a system with no difference from the first except that it has been equilibrated with water-table drawdowns consistent with historical groundwater extraction.
The natural condition is created by iteratively simulating water year 2003 (from 1 October 2002 through 30 September 2003) using hourly atmospheric forcing until the terrestrial hydrology is in approximate dynamic equilibrium. This initial condition for the coupled PF.WRF simulations was then extracted from spinup simulation results corresponding to the 0000 UTC 14 July 2003 start time.
The pumped condition was created by imposing water-table drawdown on the natural subsurface pressure field in the Central Valley portion of the domain. This was accomplished by replacing the natural conditions with a hydrostatic pressure field derived from the water level elevation simulated by the Central Valley Hydrologic Model (CVHM; Faunt 2009) for the unconfined aquifer system for the summer of 2003. The effect of this perturbation was allowed to propagate and equilibrate with the rest of the modeled domain through subsequent iterative yearly simulations using the same hourly atmospheric forcing used in the natural condition preparation. Comparison of the resulting water-table configuration to that produced by the CVHM showed the equilibration procedure did not substantively alter the agreement between the two. Water-table depths for the natural and pumped conditions are shown in Fig. 3, with contours highlighting the regions of drawdown within the Central Valley. Notably, a comparison of the plots in Fig. 3 shows the disappearance of much of the near-surface groundwater (dark blue zones) along the Central Valley axis and alluvial valleys under the modern pumped condition.
The soil moisture, defined here as the moisture in the top 2 m (four layers) in the model, was controlled by the two water-table configurations. In the natural (no pumping) condition, the soil moisture distribution mimicked the water-table depth: high soil moisture contents along the riparian areas of the San Joaquin River and main tributaries (Fig. 3), with drier soils near the Central Valley edge and to the south along the natural aridity gradient. Under the pumped water-table condition, soil moisture tended to be drier across the Central Valley floor in general, with the most intense drying occurring in the riparian zones where the natural water table was very near the surface.
d. Suite of three PBL schemes
Boundary layer processes in mesoscale models like WRF are incorporated through a parameterization that attempts to account for the subgrid-scale turbulence that cannot otherwise be explicitly simulated. A range of parameterizations is available to the modeler and each encapsulates a particular conceptual representation of the boundary layer. An ensemble consisting of simulations run with different boundary layer schemes thus permits a measure of conceptual uncertainty in the modeling result. To that end, we performed the downscaling and coupled PF.WRF simulations over the San Joaquin basin using three common boundary layer schemes: the Medium Range Forecast (MRF), the Yonsei University (YSU), and the Mellor–Yamada–Janjić (MYJ) schemes. We performed nested offline WRF simulations for each scheme choice to ensure that all initial and boundary conditions used in the coupled PF.WRF simulations were consistent with the scheme used.
These schemes are well established, commonly used, and represent different turbulence closure and local/nonlocal schemes. Specifically, the MRF (Hong and Pan 1996) and YSU (Hong et al. 2006) schemes are both first-order, nonlocal formulations. In contrast, the MYJ scheme (Janjić 1994, 1990) includes a 1.5-order closure scheme in a local formulation. The difference in underlying formulations should result in variation in simulated boundary layer characteristics. In general, the nonlocal schemes would be expected to produce deeper mixed layers and enhanced free atmosphere entrainment compared to local schemes as they capture the effects of larger eddies (Cohen et al. 2015; Stull 1991).
Small variations in initial condition can result in divergent simulated states as a result of the chaotic nature of the atmospheric system (Lorenz 1963). Given our imperfect knowledge of the initial state of the atmosphere, performing simulations with an ensemble of plausible initial conditions allows a means to assess the sensitivity of a particular result in the context of this uncertainty. For this study, we adopt the method of shifted initializations (Walser et al. 2004) to construct an ensemble of six perturbed initial conditions used as inputs to the simulation of the PF.WRF coupled domain. This method creates a set of perturbed initial conditions by 1) first initializing and running the model with initial conditions at increasing times previous to the desired start, 2) taking the difference of a given variable field and for each ensemble member from the ensemble mean at the desired start time, and 3) scaling that difference by a user-defined value and adding (or subtracting) that difference back to the selected ensemble of variable fields (Walser et al. 2004; Walser and Schär 2004). For this study we used six hourly shifts in initialization time; a scaling factor of three; and modified horizontal wind, humidity, pressure, and temperature variables. We limit the ensemble variation to the initial condition perturbations: we use identical lateral boundary conditions and the YSU boundary layer scheme in each of the six ensemble members. The same six perturbed initial conditions are used to initialize the coupled PF.WRF domain for both the predevelopment and pumped water-table scenarios.
3. Results and discussion
a. Surface energy fluxes
Changes in soil moisture are tied to shifts in the partitioning of energy fluxes at the land surface. For example, a decrease in soil moisture reduces latent heat flux and increases sensible heat flux. Given that water exists in a connected continuum in the porous subsurface, water-table depth may correlate with energy flux partitioning by way of groundwater contributions to soil moisture across a landscape. Previous studies have demonstrated an inverse relationship between water-table depth and latent heat flux (LH), with the magnitude of flux most sensitive to changes in water-table depths between 1 and 5 m (Kollet and Maxwell 2008; Rihani et al. 2010; Szilagyi et al. 2013). Such a change in LH implies a counterbalancing (at least partially) shift in sensible heat (SH) to maintain a closed energy budget. The shift in energy partitioning can be summarized by examining changes in evaporative fraction (EF), the fraction of incoming radiation that goes to latent heat flux (Gentine et al. 2011):
where λ is latent heat of vaporization, E is evaporation rate, Rnet is net radiation, LH is latent heat flux (W m−2), SH is sensible heat flux (W m−2), and B is the Bowen ratio (SH/LH). This formulation assumes that ground heat fluxes are sufficiently small that .
As a first step in understanding how water-table changes affect the system, we compute the EF for every model cell within the Central Valley portion of the domain at 30-min intervals for the 2-week simulation period. It is not uncommon to observe EF values greater than 1, particularly during deviations from daytime, clear sky conditions (Gentine et al. 2011; Peng et al. 2013). To improve clarity of the difference between the pumped and natural water-table conditions, however, we discard EF values outside the range of 0–1. The time series plots in Fig. 4 show the temporal variation of the EF over the Central Valley for the two water-table conditions and each of the three PBL schemes used. The solid lines show the spatial mean and the shaded regions represent plus/minus one standard deviation from the mean. Qualitatively, the choice to filter values to the 0–1 range does not affect the mean daytime values or the difference between the pumped and natural water-table conditions, but it does substantially reduce the standard deviation about the mean. The plots demonstrate a daily pattern in EF evolution: near zero or negative values in the evening and a positive “U shape” during daylight hours. The lowering of the water table (increasing water-table depth) leads to a decrease in the evaporative fraction during midday, with midday minima reduced from approximately 0.5 under the natural water table to 0.2–0.4 under the pumped water-table condition.
The spatial distribution of reduction in EF follows the pattern of water-table decline across the Central Valley landscape (Fig. 5). In particular, reduction in average EF is more pronounced in locations where the natural water table was nearer the land surface, like the riparian corridor along the San Joaquin River and Fresno Slough, in the southern portion of the domain. Notably, stream and river channels that maintain surface flows and/or groundwater connection under the pumped water-table condition show no change in average EF, highlighting the importance of including lateral surface and subsurface flow when studying changes across a heterogeneous landscape.
b. Boundary layer height
The simulated boundary layer follows a diurnal cycle with the amplitude and timing dependent on the location and general atmospheric conditions. The plots in Fig. 6 show a time series of the spatial mean and standard deviation about the mean of boundary layer height for the natural and pumped water-table conditions. The observed shift in EF (section 3a) suggests that the lower boundary layer receives more turbulence-generating sensible heat fluxes under the pumped water-table condition, helping to promote development of a higher boundary layer. This connection is born out in the time series of PBL heights averaged over the Central Valley: while evening and early morning PBL heights tend to match in the pumped and nonpumped scenarios, the midday average peak PBL under the lowered water-table (pump) condition is consistently greater than that under the higher, natural water-table condition. The plots in Fig. 6 also demonstrate the variability among the three PBL schemes used. The MYJ scheme produces the thinnest boundary layer but matches the diurnal pattern of the YSU while the MRF and YSU schemes produce comparable peak PBL heights while following different daily evolution patterns.
Diagnosis of the planetary boundary layer height in WRF depends on the boundary layer scheme in use. For example, the nonlocal YSU scheme defines the top of the boundary layer as the lowest elevation at which the bulk Richardson number Rib exceeds some critical value (Hong 2010; Hong et al. 2006). In contrast, the boundary layer height is defined in the MYJ scheme as the lowest elevation at which the turbulent kinetic energy falls below a specified threshold (0.1 m2 s−2; Janjić 2002). This difference in diagnosis method has been shown to increase the variance of simulated PBL height in ensembles composed of different PBL schemes (Xie et al. 2012). Acknowledging that this phenomenon may be the cause of greater between-scheme difference than a more unified measure of mixed layer depth would yield, we still present the scheme-specific diagnosed PBL height as a valuable metric of atmospheric response to change in water-table configuration.
A comparison of time-averaged peak PBL height differences (; Fig. 7) shows the change in PBL height corresponds closely with the change in EF shown in Fig. 5. Again, we see the largest increases in PBL height in locations along the Central Valley axis where water tables declined from a shallow natural depth. The plots for each of the different schemes demonstrate differences in the simulation of the PBL: the MRF scheme seems to produce a change in PBL that is most sensitive to heterogeneities in the land surface, the YSU scheme yields a PBL difference that is more diffuse but similar in magnitude to that produced in the MRF scheme, and the MYJ scheme gives the smallest change in PBL height while matching the qualitative spatial distribution of PBL change seen in the other two schemes.
The spatial plots also reveal that changes to the boundary layer extend into the Sierra Nevada to the east of the Central Valley. Changes in PBL height appear for each of the PBL schemes and are characterized by a spatial pattern of mixed PBL height increases and decreases. The apparent lack of a coherent spatial structure in these changes reflects the impact of the complex terrain on surface heating, winds, and the resulting boundary layer development. Thus, while attributing a direct mechanistic link between water-table changes in the valley to precise changes in PBL heights in the mountains may be unjustified, these results suggest the perturbation of important boundary layer properties may extend beyond the geographic location over which the perturbing subsurface change occurs.
c. Variation from PBL scheme and initial condition perturbations
The differences in simulated PBL height resulting from choice of scheme and water-table condition are shown in the time series of spatially averaged daily peak [taken as the PBL height at 1400 local time (LT)] PBL height in Fig. 8. The heavy lines show the mean PBL height response within each water-table condition group. The differences in mean PBL scheme group values indicate an average increase of 200–400 m in the peak PBL height as a result of the decline of the regional water table. The plot also shows considerable variation in the simulated PBL for the pumped water-table condition that depends on the PBL scheme, with the mean peak PBL for the MYJ scheme roughly matching that of the natural water-table scheme ensemble mean. We test the significance of the differences in PBL height due to the water-table configuration, choice of scheme, and the interaction effect of these two factors, using a repeated measures analysis of variance (ANOVA). The daily spatial mean peak PBL heights for the Central Valley portion of the domain (n = 15) were used as the repeated measures within each PBL and water-table group. The results of the ANOVA (Table 2) show that both the water-table depth condition and the choice of PBL scheme have a significant effect on the resulting peak PBL height at the p < 0.001 level. A small interaction effect is also observed, with a significance level of roughly 0.05. This means that the choice of PBL scheme has some effect on the magnitude of the difference in average peak PBL height resulting from different water-table depths, although the strength of this effect appears to be less than the major effects of PBL scheme or water-table condition alone.
The differences in simulated daily peak average PBL heights resulting from initial condition perturbations are shown in Fig. 9. In contrast to the effects of PBL scheme (Fig. 8), the variation in peak PBL height during the simulation due to initial condition perturbation are minimal. An analysis of variance shows the lowering of the water table has a significant impact (<0.001) on average PBL height while the effect of initial condition perturbation is insignificant. Although the initial condition perturbations had almost no effect on the PBL, they resulted in variation in other parts of the system. For example, domain-averaged accumulated precipitation varied between 2.59 and 2.84 mm in the natural water-table condition and between 2.62 and 2.85 mm in the pumped water-table condition. The coefficient of variation is approximately 3% in both cases and is consistent with results from Walser and Schär (2004) for cases of moderate convection and strong orographic control.
d. Comparison to observations
Comparing the results of simulations in this study against observations is complicated by two factors. First, the hydrology in the modeled system reflects either a natural groundwater system or one affected by water-table drawdown absent accompanying irrigation at the land surface. The atmospheric measurement record, insofar as it exists, captures the effects of both water-table drawdown and widespread irrigation. Second, measurements that overlap with the spatial extent of the domain and the simulation period are relatively scarce. A station featuring collocated radar wind boundary layer profiler (915 MHz) and radio acoustic sounding system (RASS) at Chowchilla [station identification: CCL, NOAA/ESRL/Physical Sciences Division (PSD); https://www.esrl.noaa.gov/psd/data/obs/datadisplay/] provides the only measure of boundary layer and near-surface dynamics over the July 2003 simulation time period. Hydrologic inconsistencies notwithstanding, a qualitative comparison between profiler estimates and the simulated boundary layer profile (Fig. 10) provide some confidence that the coupled model correctly captures some aspects of lower atmosphere dynamics.
The comparison of measured and simulated (YSU PBL scheme and natural water-table condition) virtual temperature in the boundary layer at three times in the diurnal cycle (Figs. 10a–c) shows that the model and observations are generally consistent within the range of observed and simulated variability. Given the impact of irrigation is embedded in the boundary layer measurements and that such activities are not simulated, an exact match would not be expected. Figure 10d shows a comparison of simulated boundary layer height (YSU–no pump, solid line; YSU-pump, dashed line) with the signal-to-noise ratio obtained from the radar profiler, a measurement that can provide an inexact, qualitative assessment of boundary layer height (Knupp et al. 2006; Heinselman et al. 2009). The comparison shows the model matches the diurnal cycle of PBL development fairly well, but tends to overestimate boundary layer height. This discrepancy would be expected given the boundary layer sensed by the radar profiler would have been affected by widespread irrigation (with the potential to dampen PBL growth) while the simulated scenarios reflect the natural seasonal (drier) soil moisture.
e. Water table–PBL relationship
The complexities of the Sierra Nevada and Central Valley landscapes represented in the PF.WRF Model introduce variation that helps define the characteristics of the water table–PBL relationship. First, the terrestrial hydrologic system simulated by ParFlow yields a spatially varying water-table configuration. Lateral subsurface gradients drive groundwater toward the riparian discharge zones along the axis of the valley while higher elevations (valley edges and the mountain block) serve as deep water-table recharge zones. Second, water-table drawdowns, as represented by the CVHM simulation results (Faunt 2009), are not uniform in space but rather represent the combined result of hydrology, climate, soils, and human decision-making. The combination of these two factors yields a system characterized by a range of hydrologic change relative to the natural state but bounded by the dynamics of the system. This allows an examination of the sensitivity of the simulated PBL response to water-table changes using each cell as a sample from the system distribution.
The plots in Fig. 11 show the relationship between water-table depth and peak (1400 LT) PBL height, averaged over the 2-week simulation period, for the natural [Fig. 11 (top) and blue points in Fig. 11 (bottom)] and pumped [Fig. 11 (middle) and red points in Fig. 11 (bottom)] water-table conditions and for each of the PBL schemes Fig. 11 (left, center, and right). Locations where the water table intercepted the surface (wetland and stream channels) were removed from the plot so that only positive nonzero values for water-table depth remained. Figures 11 (top and middle) show a two-dimensional histogram (or alternately a “scatter density” plot, with color intensity indicating higher frequency) plotted with the water table (the x axis) on a logarithmic scale. Figure 11 (bottom) overlays the data points from Figs. 11 (top and middle) and plots the water table on a linear scale.
The plots in Fig. 11 indicate a variable relationship between water-table depth and PBL height. In general, the peak PBL height increases with increasing water-table depth, although much of this change occurs where water tables are more than several meters deep. From a visual inspection of the PBL height–water-table depth plots, the relationship can be categorized into three regimes. The first regime exists over very shallow water tables and is characterized by generally lower but variable PBL heights. Based on the results presented here, we define this regime to exist over the range of water-table depths from near zero to 7 m. This range is somewhat arbitrary but is chosen here as being consistent with a visual “break” in the scatter density plots. Interestingly, there appears to be little overall correlation between water-table depth and boundary layer thickness over this water-table depth range, as evidenced by the nearly horizontal clustering of points in the 2D semilog histograms and the large spread of points in the linear scatterplots. The second regime is characterized as a zone of approximately log-linear correlation between water-table depth and peak PBL height and extends from a water-table depth of 7 m to approximately 150 m. The envelope of the scatter points tends to condense over this range toward a common positive trend. The third regime exists where peak PBL heights reach an asymptotic limit and effectively decouple from boundary layer processes, at water-table depths of approximately 150 m. This regime is the least thoroughly sampled of the three (given the relative rarity of such deep water tables in this domain) and thus is less discernible in the scatter density plots but evident in the linear scatterplots.
The variation in behavior between the MYJ scheme and the YSU and MRF schemes presents challenges to generalizations across schemes. The plots for the YSU and MRF schemes indicate the PBL is more sensitive in the intermediate water-table depth regime compared to the MYJ scheme. Given that the shift in EF is similar in all schemes, this difference in boundary layer sensitivity highlights the effect of a local PBL scheme such as the MYJ.
The increase in average peak PBL height with increasing water-table depths beyond 10 m is not fully explained by the results of previous studies that suggest relative insensitivity of surface fluxes to local water-table depths exceeding 7–10 m (Kollet and Maxwell 2008; Maxwell and Kollet 2008; Szilagyi et al. 2013). However, the separation of points and shift in log-linear slope between the pumped and natural water-table scenarios (Fig. 11) over this water-table depth range (7–150 m) in this study suggests that the regional configuration of the water table determines the intensity of the water table–boundary layer relationship. Put another way, the frequency distribution of water-table depths across the simulated Central Valley region appears to affect how the boundary layer responds to local land surface–subsurface connections. As an example, consider that a comparison of the points at a given water-table depth, say at 20 m, shows that the average peak PBL height value will tend to be greater for the pumped water-table condition than for the natural water table.
The PBL response to subsurface moisture, then, may be dependent on a combination of local and regional conditions. A wetter landscape, as measured by the higher proportion of locations with a shallow water table, may suppress PBL development even in locations with deep water tables (and locally dry soils) as a result of mixing of cool, moist near-surface air within the boundary layer or the production of secondary circulations that moderate boundary layer development, as in Shrestha et al. (2014). Conversely, as regional groundwater levels decline and the landscape becomes drier, net radiation will be partitioned toward sensible heat flux over a larger portion of the land surface, thus reducing the spatial extent and/or intensity of any inhibiting effect on PBL development over comparably dry soils. Furthermore, the lowering of the water table in the Central Valley tends to make soil moisture more uniform (drier) over larger areas, leaving fewer patches of adjacent wet and dry soil that can generate boundary layer inhibiting circulations.
The scatterplots reveal a complex dependence of peak PBL height on local and regional water-table depth. We can examine this relationship further by calculating the change in peak PBL height between the pumped and natural water-table scenarios relative to the corresponding change in water-table depths. Plotting this variable against the pumped water-table depth summarizes the sensitivity of the boundary layer as a function of the water-table state (Fig. 12). Despite considerable scatter in the values, the plotted points define an envelope of responses that indicates the PBL is most sensitive to change in water-table depth when that change leads to a water table that is less than 10 m deep. This lower sensitivity limit is consistent with the idea that vertical redistribution of water from the phreatic surface becomes negligible once the water table reaches some depth threshold of 5–10 m, the so-called “critical zone” (Kollet and Maxwell 2008; Maxwell and Kollet 2008). Notably in contrast to this critical zone concept, however, is the absence of an upper water-table limit, that is, a reduction in sensitivity at shallow water-table depths. Rather, the PBL appears to be the most sensitive when water-table depths change within 5 m of the land surface. This may be a consequence of the midsummer simulation period used in this study. The midday evaporative demand of the boundary layer is always sufficiently high that evaporation is rarely energy limited, at least on average across the simulated Central Valley landscape. This ensures that a change in shallow water-table depth sufficient to change the ET-supporting vertical flux will necessarily reduce ET and shift surface energy partitioning toward PBL-enhancing sensible heat flux. Whether the absence of an upper limit on boundary layer sensitivity to water-table changes exists for different seasons or locations remains a question for future study.
f. Dynamic boundary layer response
As a final analysis of boundary layer response to changes resulting from a regionally lowered water table, we briefly examine two measures of boundary layer dynamics: vertical wind and turbulent kinetic energy (TKE). Figure 13 shows temporally averaged vertical cross sections of vertical wind, TKE, and soil saturation along an east–west profile through the middle of the model domain (model row 120). Note the vertical axis by model layer and model component is not to scale. Both plots show that reductions in soil moisture (driven by a lowered water table) propagate to boundary layer as increased turbulence and changes in vertical wind patterns. In the case of TKE, the effect is uniform: across the Central Valley floor the TKE increases with a lowered water table and appears to increase more in association with intense clusters of soil moisture reduction. In contrast, average vertical wind differences are not spatially uniform but rather reflect the patterns of soil moisture difference. Locations with more intense reductions in soil saturation tend to be associated with an increase in mean vertical wind speed, either more strong upward wind or weaker downward wind. Likewise, less intense reductions in soil moisture are associated with stronger downward wind or weaker upward wind. Because spatial heterogeneity in soil moisture driven by groundwater extraction has been shown to instigate local circulations in the boundary layer (Shrestha et al. 2014), the results shown here further suggest a mechanism by which local variability may moderate (or enhance) a regionally dominant change in boundary layer properties.
We simulated the connected terrestrial and atmospheric hydrologic system of the San Joaquin River basin in central California using a coupled ParFlow–WRF Model. Using this simulation framework, we tested the impact of a regional water-table decline, consistent with historical groundwater extraction in the Central Valley, on the planetary boundary layer for a 2-week period in July 2003. Spatial mean boundary layer heights over the Central Valley increase with the water-table decline. This effect is statistically significant in results simulated using three common PBL schemes and six initial condition perturbations, although the choice of PBL scheme appears to affect the magnitude of the mean PBL height difference between the natural and pumped water-table conditions. The relationship between peak daily PBL and water-table depths at individual locations is complex; results suggest the correlation depends on the water-table depth regime, the regional water-table configuration (natural vs historic drawdown), and choice of PBL scheme.
Scatterplots reveal that the PBL is highly variable over extremely shallow water tables, although peak PBL heights tend to be lowest for water-table depths less than ~7 m. Over locations of deep water tables (>150 m) peak PBL height tends asymptotically toward a maximum value, indicating a decoupling between groundwater and boundary layer processes. In zones of intermediate water-table depths, that is, 7–150 m, the relationship between PBL height and water-table depth is approximately log linear, suggesting a weak coupling between the saturated subsurface and the atmosphere. Each of the schemes tested yields qualitatively similar relationships, but the local MYJ scheme seems to dampen the intensity of the PBL–water table relationship compared to the nonlocal schemes. In general, lowering the regional water table, as we did here by imposing historic drawdowns on an approximation of the natural system, increases the overall slope of the PBL–water table relationship. In other words, one would expect higher peak PBL heights (compared to a higher water-table configuration) for cases of regional water-table drawdown even in locations with equivalent water-table depths. Furthermore, the sensitivity of the PBL–water table relationship is highest for water-table changes within the top 10 m, despite this region being comparatively less sensitive to water-table depth for a given regional water-table configuration.
This study demonstrates the connections between groundwater and the atmosphere for a real-world system. However, the range of atmospheric and terrestrial conditions tested limits a broader interpretation of this study’s conclusions. Variation in simulation factors such as time of year, location, and domain extent would be valuable in testing the robustness of these findings. Moreover, including the asynchronous effects of irrigation on soil moisture in coupled groundwater–atmospheric simulations would provide insight into the relative role of moisture distribution throughout the subsurface on land–atmosphere interactions and feedbacks. Finally, the 1-km lateral and variable vertical ParFlow model discretization used in this study represent one of many possible discretizations that could be used. More finely resolved terrestrial hydrology may result in more heterogeneity in soil moisture and surface fluxes that can, in turn, affect boundary layer development (Shrestha et al. 2015). Finer vertical discretization in the subsurface would permit a more detailed simulation of water-table dynamics during transient drawdown periods (e.g., during an irrigation season). Such refinements in lateral and vertical subsurface resolution and their impact on groundwater–atmosphere interactions are important considerations that warrant further study.
Despite these limitations, the results of this study have implications for understanding the propagation of impacts associated with anthropogenic groundwater extraction through the connected hydrosphere. First, this study’s conclusions build upon those in previous work and provide further evidence in support of including a physically consistent representation of deeper subsurface water in coupled land–atmosphere simulations. That the variable relationship between peak boundary layer height and water-table depth extends on the order of 10 m below the land surface suggests that a more vertically limited representation of subsurface hydrology may fail to capture the range of system interactions. Second, water-table declines are not isolated to California’s Central Valley but are a global phenomenon, with regions of intense irrigation and extraction continuing to deplete aquifer systems (Döll et al. 2012; Scanlon et al. 2012; Wada et al. 2010, 2012). Given the spatial extent of water-table decline, groundwater effects on the planetary boundary layer may represent another anthropogenic forcing that warrants consideration for studies of regional weather and circulation.
We would like to acknowledge high-performance computing support from Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. This research was supported by funding from the National Science Foundation through its Climate Change, Water, and Society (CCWAS) Integrative Graduate Education and Research Traineeship (IGERT) program (http://ccwas.ucdavis.edu/; DGE-1069333).
Current affiliation: Technical Service Center, U.S. Bureau of Reclamation, Denver, Colorado.