2.4.1. ParFlow-CLM

ParFlow is a three-dimensional integrated hydrological model that simulates subsurface and surface water flows by solving the Richards’ equation and shallow surface water equations [38–41]. ParFlow contains a coupled land surface module, CLM, which solves the energy balance for many land surface processes. Canopy water balance, losses and additions from evapotranspiration (ET), precipitation and snowmelt are communicated with ParFlow at every timestep [33, 34, 42–44]. ParFlow-CLM has been applied to the ERW at 1km and 100m resolution; full details of the model construction and performance can be found in Foster and Maxwell [35] and Foster et al. [45]. Here, we use both input and output from the 100m resolution model run. The subsurface is discretized into five layers across the entire watershed, the top three are soil layers of depths 0.1m, 0.3m and 0.6m, while the bottom two are geological layers of 8m and 21m, with the deepest 21m representing fractured bedrock.

Pressure head output from ParFlow-CLM was spatially integrated for each sub-watershed i (Fig 1) to determine the total volume of water stored in the groundwater, V_g,i, vadose zone, V_s,i, and surface water, V_r,i, at hourly intervals. Movement into and out of the vadose zone is quantified with two fluxes from ParFlow-CLM, also at hourly intervals: infiltration into the vadose zone, Infilt_i, and evapotranspiration from the soil layers, ET_s,i. Infiltration includes snowmelt water, precipitation, and runoff from adjacent cells that enters the vadose zone, while exfiltration (i.e. negative infiltration) includes vadose zone water that exits to the surface via saturation excess. Surface water pressure head values were converted to discharge (flows) using Manning’s equation for the outlet of each sub-watershed. Discharge is a function of the representative slope values of each outlet, and the Manning n value of that cell, as parameterized by land cover type. As discussed in Foster and Maxwell [35], Manning n and hydraulic conductivity were used constrain system-wide discharge by manual calibration with stream-level observations to control the dynamics of the streams.

Soil and air temperature (T_soil and T_air, °C) were output as hourly spatial averages for each sub-watershed. Air temperatures are derived from PRISM datasets [46], interpolated to hourly resolution using phase two of North American Land Data Assimilation (NLDAS-2) forcing [47], which are available at 1/8^th-degree at hourly time steps, and were interpolated and downscaled to match the discretization of the ParFlow-CLM model. Soil temperature is solved in CLM using the heat diffusion equation and a subsurface heat flux with the Fourier law for heat conduction over the top 2 meter of the model for both soil and snow layers [48]. Stream water temperatures were calculated using the empirical relationship given in Lauerwald et al. [49]. We calculate a single groundwater temperature (T_gw, °C) for the entire ERW by averaging sub-watershed soil temperatures at each timestep.