Figure 2. The elevation (a), land cover (b), and soil (c) and aquifer (d) properties in the modeling area.
  1. Discretization and coupling depth
The subsurface is divided into 1 × 1 km2 cells in the horizontal direction with terrain following grids to accommodate the topography. Layering is customized in this study. The subsurface of 100 m thickness is discretized into 15 layers with 6 top layers of 0.1, 0.3, 0.6, 1.0, 3.0 and 5.0 m and 9 bottom layers of 10 m per layer. Thus, the number of grid cells are 25215 (41 × 41 × 15) in total. The top 4 layers (2 m) are set to be soils (Figure 2c) while the deep subsurface (98 m) has the hydraulic properties of the aquifers (Figure 2d). Six scenarios with different coupling depths (2, 10, 20, 30, 60, and 90 m) are tested to study the feedbacks between land surface and subsurface processes, especially with the long-term GW pumping. The coupling depth corresponds to the number of coupling layers of 4, 6, 7, 8, 11, and 14, respectively. Hence, the scenarios are noted as L4, L6, L7, L8, L11, and L14 respectively, where L represents the layer.
  1. Boundary conditions and simulations
The boundary conditions are all no flow boundaries except for the overland flow boundary at the land surface. Meteorological forcings from 10/01/1970 to 09/30/1971 in this study are obtained from the 1.25-degree gridded Japanese 55-year Reanalysis (JRA-55) (Harada et al., 2016; Kobayashi et al., 2015). The forcings are spatially uniform by using the data near the center of the modeling area (Kollet and Maxwell, 2008). The original JRA-55 data of 3-hour resolution are linearly interpolated into hourly resolution. The 1-year forcings are repeatedly used in the following 10-year pumping simulations.
After spin-up (refer to section 2.2.4), one more year is simulated to represent the natural state without pumping while 10 years of simulation with pumping is conducted to explore its long-term effects on GST. GW pumping is assumed only in the croplands shown in Figure 2b. Multi-year averaged annual pumping rate of about 80 million cubic meters is adopted from Condon and Maxwell (2014a). One more group of pumping simulations with tenfold pumping rate are also conducted to understand the mechanisms of the subsurface buffer. Those two groups of simulation are noted as G1 and G2 respectively. Hourly time step is used in all simulations.
  1. Initial conditions
Initial conditions are critical and specifically considered in this study due to the thickened coupling depth. In Stevens et al. (2007), a spin-up of 500 years was conducted for 500 years of simulation, and the heat propagating downward did not reach the 1000 m BBCP. In this study, for L4, which is the common set in previous studies (e.g., Maxwell and Condon, 2016), 4 years of spin-up are enough to achieve the dynamic equilibrium of the model. Nevertheless, considering the 10 years of pumping, longer spin-ups of 10, 30, 50, and 100 years (≥ 10 years pumping) were performed additionally to ensure the reliability of the simulations. Therefore, the following discussion in section 3 is based on 4 years of spin-up while the results with longer spin-up are provided in Supporting Information (Figures S1 and S2).
  1. Model evaluation
The model was evaluated in Condon and Maxwell (2014a), which demonstrated its capability to capture the interactions of water and energy processes between land surface and subsurface. Also as pointed by many previous studies based on such integrated models (e.g., Condon and Maxwell, 2014b, Maxwell and Condon, 2016), the aim here is to diagnose the subsurface buffer based on a realistic platform but not to fit history or predict future by a calibrated model. Hence, no special calibration and verification of the model were conducted in this study.
  1. Results and discussion
  2. Performance of the subsurface buffer
  3. WTD and GST variations under long-term GW pumping
If not otherwise specified, the discussion in this section is based on scenario L4 in G1. Annual average WTD in the year without pumping is 2.32 m in average and 25.08 m at the maximum (Figure 3a). The difference between annual average WTD in each pumping year and that in the year without pumping was calculated (ΔWTD) (Figure 3). Obvious increase of WTD occurs not only in the croplands with pumping but also at the uplands without pumping (Figure 3), which is consistent with Ferguson and Maxwell (2011) and Condon and Maxwell (2014a). The mean ΔWTD in the modeling area was plotted in Figure 4. The mean ΔWTD after 10 years of pumping is about 0.5 m for G1 and 10 m for G2. Hence, after 10 years of pumping, for G1, most WTD in the modeling area is still in the critical depth range (1–10 m), in which the land surface heat fluxes are sensitive to the variations of WTD (Condon and Maxwell, 2019; Ferguson and Maxwell, 2011). In contrast, for G2, most WTD has become lower than the critical depth range. In addition, the rate of increase of WTD is decreasing with time in G1 while it is almost linear in G2. The mean, maximum, and minimum change of GST (ΔGST) are shown in Figure 5. GW pumping leads to increased GST in summer (March to September) and decreased GST in winter (October to February) (Figures 5a and 5c). The mean ΔGST indicates the warming trend in average (Figures 5a and 5c), though the decrease of GST is over its increase at extreme values at some moments (Figure 5b). The mean ΔGST is in -0.5 K–0.5 K in G1, and in -1.5 K–1.5 K in G2. In addition, the variations of GST occur most obviously in the pumping area (Figure 6). These verify the hypothesis proposed in section 1 (i.e., the first objective in this study) that the subsurface may be conceptualized as a buffer on variations of GST in the Little Washita basin. GW pumping can weaken this buffer by causing hotter summer and colder winter with a warming trend in average.