Numerical Terradynamic Simulation Group 5-2013 Future humidity trends over the western United States in the CMIP 5 global climate models and variable infiltration capacity hydrological modeling system

Global climate models predict relative humidity (RH) in the western US will decrease at a rate of about 0.1– 0.6 percentage points per decade, albeit with seasonal differences (most drying in spring and summer), geographical variability (greater declines in the interior), stronger reductions for greater anthropogenic radiative forcing, and notable spread among the models. Although atmospheric moisture content increases, this is more than compensated for by higher air temperatures, leading to declining RH. Finescale hydrological simulations driven by the global model results should reproduce these trends. It is shown that the MT-CLIM meteorological algorithms used by the Variable Infiltration Capacity (VIC) hydrological model, when driven by daily Tmin, Tmax, and precipitation (a configuration used in numerous published studies), do not preserve the original global model’s humidity trends. Trends are biased positive in the interior western US, so that strong RH decreases are changed to weak decreases, and weak decreases are changed to increases. This happens because the MT-CLIM algorithms VIC incorporates infer an overly large positive trend in atmospheric moisture content in this region, likely due to an underestimate of the effect of increasing aridity on RH. The result could downplay the effects of decreasing RH on plants and wildfire. RH trends along the coast have a weak negative bias due to neglect of the ocean’s moderating influence. A numerical experiment where the values of Tdew are altered to compensate for the RH error suggests that eliminating the atmospheric moisture bias could, in and of itself, decrease runoff up to 14 % in high-altitude regions east of the Sierra Nevada and Cascades, and reduce estimated Colorado River runoff at Lees Ferry up to 4 % by the end of the century. It could also increase the probability of large fires in the northern and central US Rocky Mountains by 13 to 60 %.


Introduction
The Earth's climate is warming due to the accumulation of human-produced greenhouse gases in the atmosphere (IPCC, 2007). Over the oceans, warmer surface temperatures will likely lead to increased evaporation and therefore greater specific humidity, but an approximately constant relative humidity (RH); the greater concentration of water vapor will in turn warm the surface further, since water vapor is a potent greenhouse gas (GHG) (e.g., IPCC, 2007;Dessler and Sherwood, 2009). This feedback mechanism depends on warmer planetary temperatures leading to increased atmospheric water vapor concentrations, which has already been identified in satellite observations (Santer et al., 2007) and surface humidity measurements (Willett et al., 2007).
It seems unlikely that RH will remain constant in locations far from large open bodies of water and where annual evaporative demand considerably exceeds annual precipitation, such as the arid regions of the western US. In such areas, warmer surface temperatures, along with limited moisture availability, may lead to lower relative humidities in the future than are experienced today. This could have an effect on hydrological and ecological processes that are sensitive to humidity, such as evapotranspiration (Friend, 1995), runoff, wildfire (Brown et al., 2004), and plant growth (Leuschner, 2002). Irrigation can increase surface relative humidity locally (Kueppers et al., 2007), but we do not include the effects of irrigation in this work.
A large number of studies have examined changes in the hydrology of the western US and related processes using the Variable Infiltration Capacity (VIC) hydrological model (Liang et al., 1994). Hydrological processes are strongly influenced by topography, so VIC is often run with a resolution of 1/8 × 1/8 • (longitude by latitude) or finer. Even with the increasing resolution of global climate models, which can approach 1/2 by 1/2 • , VIC runs are likely to be used for hydrological applications for a considerable time yet. VIC can also easily be run over a limited area, increasing efficiency further.
VIC is a well-regarded model ( (Liang et al., 1994) has been cited nearly 600 times according to the Thomson Reuters Web of Knowledge) that does a good job of simulating historical changes in streamflow when driven with historical meteorology (e.g. Abdulla and Lettenmaier, 1997;Hamlet and Lettenmaier, 1999b;Arnell, 1999;Nijssen et al., 2001b;Wang et al., 2008;Niu and Chen, 2009), and has been used extensively both in the region and globally. For example, Hamlet and Lettenmaier (1999a) used VIC for forecasting Columbia River streamflow based on the patterns of climate variability; Nijssen et al. (2001a) examined the sensitivity of global river flow to climate change; and VanRheenen et al. (2004) used VIC and downscaled global climate model simulations to study changes in flow in the Sacramento and San Joaquin basins. Christensen et al. (2004) and Christensen and Lettenmaier (2007) did similar analyses for the Colorado River basin, while Hayhoe et al. (2004) examined future changes in California water resources. Westerling et al. (2006), Westerling and Bryant (2008), and Westerling et al. (2011) used VIC simulations to help understand historical trends and future projections in wildfire activity, and Hamlet and Lettenmaier (2007) did a similar analysis for runoff, evapotranspiration, and soil moisture. Maurer (2007) examined runoff over the Sierra Nevada; Barnett et al. (2008) and Pierce et al. (2008) used VIC to study changes in runoff, streamflow, and snowpack over the western US; and Adam et al. (2009) examined impacts on snowpack globally. VIC has been used for a wide variety of purposes in many locations.
The studies enumerated above are driven by daily minimum temperature, maximum temperature, and precipitation (T min , T max , and P ). Often, these fields are derived from monthly global climate model output that has been resampled to a daily timescale using a statistical downscaling technique , although statistically downscaled daily global model data has been used as well . Wind speed is another required input field, and is sometimes specified using climatological values (e.g., from Maurer (2002)), since statistical downscaling procedures have typically not yet been tested with wind speed.
The hydrological simulation needs additional meteorological variables, such as humidity and incoming shortwave radiation. Since these are relatively sparsely observed, VIC incorporates algorithms taken from the MT-CLIM version 4.2 package (Hungerford et al., 1989;Thornton and Running, 1999;Thornton et al., 2000; see also Maurer, 2002) to calculate these fields from T min , T max , and P . A global evaluation of the MT-CLIM algorithms can be found in Bohn et al. (2013). Henceforth we will call this combination of the VIC hydrological model and the MT-CLIM meteorological algorithms the VIC modeling system (VMS). We introduce this terminology to emphasize the distinction between VIC (the hydrological component of the system) and the entire modeling package that enables hydrological simulations to be computed from a given set of meteorological fields. Note, however, that the references given above do not draw this distinction, and commonly identify hydrologic simulations derived using the MT-CLIM algorithms and VIC as "VIC" simulations.
Checking VMS's ability to simulate humidity changes is important given the large number of published climate change studies that have analyzed streamflow or runoff using this approach. Results from VIC have informed our understanding of the effects of climate change on runoff, streamflow, wildfire, and snowpack in the western US and around the world, and could potentially influence key resource, policy, or adaptation decisions.
VIC uses humidity when calculating the incident solar radiation (along with solar angle, elevation, slope, etc.) and evaporation from bare soil. Additionally, when the relative humidity drops, modeled stomatal resistance to transpiration in plants increases; this effect is observed in nature and included in VIC's parameterizations. Under high relative humidity conditions the model includes the possibility of condensation and dew formation on the ground, which increases the soil moisture, but this is a small term and in most places represents a negligible contribution to the soil moisture budget.
The MT-CLIM meteorological algorithms as incorporated in VMS assume that the dew point temperature, T dew , is closely related to nighttime minimum temperature. Although T dew is not an ideal proxy for atmospheric moisture content because it depends in part on atmospheric pressure (and thus elevation) and is nonlinearly related to moisture content, it nevertheless is a good moisture indicator and has the advantage that it is more commonly observed at stations than other humidity variables (Robinson, 1998;Brown and DeGaetano, 2009). In coming decades T min will increase at least partly due to greater concentrations of CO 2 and other non-water vapor GHGs, yet VMS equates an increase in T min with an increase in humidity. There is a danger that VMS may overpredict T dew increases as a result.
In humid regions the assumption that T min equals T dew is supported by observations, especially in clear and calm conditions (Brown and DeGaetano, 2009). However in arid regions (such as much of the western US), the dew point is often lower than T min . Accordingly, the meteorological algorithms in VMS incorporate a correction to the T dew calculation suggested by Kimball et al. (1997)  , which takes into account the local aridity. The Kimball-97 T dew correction is empirical in nature, being based on a polynomial fit to observations. It might therefore become a progressively worse approximation as anthropogenically produced greenhouse gases alter the radiation balance of the planet.
A number of recent studies have focused on current and future humidity trends. Much of this work has explored the dynamics and thermodynamics of humidity changes in the atmosphere (e.g., Sherwood et al., 2010;Schneider et al., 2010;Seager et al., 2010;Wright et al., 2010). Other, observationally based work has examined changes in surface relative humidity. Gaffen and Ross (1999) found that relative humidity increased in the northern interior US from 1961-1995, but slightly (non-significantly) decreased during summer and autumn in the southwest. Dai (2006) analyzed global weather reports from 1976 to 2004 and found RH increases in the northern part of the western US on into Canada, and decreases in the southern part, similar to Gaffen and Ross (1999). On the other hand, Willett et al. (2007) found RH declines in this entire region over the period 1973-2003, although the trends were small and statistical significance was not indicated. Vincent et al. (2007), analyzing station data over the period 1953-2005, found generally negative but not statistically significant RH changes in the interior western part of Canada, roughly in accord with Willett et al. (2007). Robinson (2000) found more complicated patterns in T dew trends that depended on season and time of day.
The purpose of this work is to show projections of future humidity changes over the western US from the new CMIP5 (Taylor et al., 2012) set of global climate models, how the trends are misestimated by VMS, and the implications this has for future runoff changes. We start with a description and evaluation of the performance of VMS's humidity simulation given historical observations of T min , T max , and P (Sect. 2). We then show future changes in surface RH over the western US projected by global climate models (Sect. 3.1), and evaluate how well these changes are captured by VMS when it is driven by T min , T max , and P from the global model simulations (Sect. 3.2). Global climate model fields are often bias corrected and/or downscaled before being used to drive a hydrological model, so we examine the effect of those processes on the RH trends as well. We additionally run a test simulation with VMS to see what effect the humidity errors might have on runoff (Sect. 4). A summary discussion and conclusions are given in Sect. 5.

Estimating humidity from temperature and precipitation
Our primary interest is in future RH trends, which are determined by both air temperature and atmospheric moisture content trends. VMS is supplied with air temperatures, but must estimate atmospheric moisture content, so any errors in the RH field that can be attributed to VMS must arise from errors in atmospheric moisture content. VMS' atmospheric moisture content algorithms estimate T dew , and T dew is available from meteorological stations, so it is a natural choice for comparing to observations. The following terminology will be used to describe the T dew calculation. The input variables that the meteorological algorithms in VMS use are daily T min , T max , and P , which need to be supplied from an external source such as observations or a global climate model. The Kimball-97 algorithm is a function of three parameters, which are EF , DTR, and T min . In Kimball-97's original formulation, EF is the nondimensional ratio of daily potential evapotranspiration (PET) to annual precipitation, and so is a measure of aridity (higher values mean greater aridity). DTR is the diurnal temperature range, T max − T min . Finally, T min is both an input variable and a parameter.
VMS converts the three input variables into the three parameters required by Kimball-97 using the MT-CLIM algorithms (Hungerford et al., 1989;Thornton and Running, 1999;Thornton et al., 2000). The conversion of T min and T max into DTR is trivial, and of course T min is used without conversion. EF is computed by first estimating PET using the Priestley-Taylor equation (Priestley and Taylor, 1972), then dividing by annual P . Of course, the Priestley-Taylor PET estimates themselves use the simplified radiation estimates derived from T min and T max , so estimated PET likely has errors that subsequently contribute to errors in EF . The version of the MT-CLIM algorithms used in VMS employ a modified version of Kimball-97's T dew correction that computes EF using the annualized precipitation over the past N days rather than the annual precipitation (i.e., using [P summed over the last N days] × [days per year/N]). Additionally, the annualized precipitation is not allowed to drop below a minimum value, P min . This introduces two new parameters: N (default: 90 days) and P min (default: 8 cm). Below it will be shown that this modification has a detrimental effect on the simulation of T dew in Mediterranean coastal locations such as central and southern California. The VMS input files allow one to specify each grid cell's annual average P , but this specified value is not used in the humidity calculation in versions of VMS released since late 2004 (we do not have access to earlier versions).

Kimball et al.'s T dew parameterization
The Kimball-  Eq. (4): where temperatures are in Kelvin, DTR is that day's diurnal temperature range and EF is the dimensionless ratio of that day's potential evapotranspiration to annual precipitation (in VMS, to 90-day annualized precipitation subject to a minimum of 8 cm). Figure 1 shows observed daily values of the left-hand side of Eq. (1) (dots) as a function of EF for 6 selected stations, 3 along the west coast and 3 in the arid continental interior. Also shown (line) is the right hand side of Eq. (1), i.e., the Kimball-97 suggested correction to T dew as a function of aridity. The detrimental nature of the suggested correction in locations influenced by the marine environment is obvious. In places such as San Diego, Long Beach, and Oakland, even during days of high evaporative demand (large EF ), the essentially infinite moisture source of the ocean keeps T dew /T min near 1. This is expected, but not captured in the correction. MT-CLIM exhibits errors in downward shortwave radition in the coastal environment as well (Bohn et al., 2013).
Perhaps more surprising is that even in the dry interior, the suggested correction, which is cubic in EF , is inferior to a simpler linear fit. The r 2 values for the Kimball-97 fits in Fig. 1 are 0.04, 0.01, and 0.18 for El Paso, Tucson, and Las Vegas, respectively, while the r 2 values for a linear fit are 0.24, 0.30, and 0.40. This is seen at many locations, not just those shown in Fig. 1. Also, Kimball-97's parameterization consistently overestimates T dew at EF = 0, which may arise from a requirement that T dew equals T min when EF equals 0. Such a requirement is not supported by the data shown in Fig. 1, or more generally in data from the other western US stations used here. However, it should be noted that data from humid east coast locations was not examined, and could display different characteristics. Figure 1 suggests that the T dew correction could be improved straightforwardly. The correction could be recast to a linear fit in EF . The endpoint at EF = 0 could be allowed to take on a value other than exactly 1.0, consistent with the observations. And instead of the cubic coefficient, which is not required if a linear correction is used instead, a new parameter could be introduced that would reflect whether or not the location was influenced by the marine environment. For example, the distance to the oceans, or even a simple binary flag for coastal vs. continental environment. Such improvements in the algorithm would be particularly valuable when using VMS in coastal regions.

Performance of the humidity algorithm in the western US
In this section the ability of VMS to estimate humidity (specifically, T dew ) across the western US will be evaluated using daily observations of T min , T max , P , and T dew from 74 meteorological stations across the western US. The first three observed variables will be used to drive VMS, and VMS's estimated T dew will then be compared to the observed value. Seventy-four global summary of day (GSOD) stations across the western US were selected for this analysis (Fig. 2), and their daily data downloaded from ftp://ftp.ncdc.noaa. gov/pub/data/gsod/. The period included was 1975 through 2009. The stations are a superset of the hourly data stations used in Kimball-97, but concentrated in the western part of the continent. (Note that in Kimball-97, only one year of hourly data was used for each station, with different years for different stations.) The stations reflect population patterns rather than topography and climate patterns, and so do not sample interior snow-producing mountain ranges or high elevations in proportion to their importance to generating runoff. The GSOD values are based on the GSOD reporting period of a UTC based day, while calculation of daily T min from the hourly values is based on a local-time day. To examine if this reporting-period discrepancy affected our analysis, we compared the T dew and T min time series calculated from the hourly and daily data at 8 stations (KSAN, KTUS, KRDM, KLGB, KOAK, KELP, KLAS, selected to include both inland and coastal locations) over the period 1980-2005, and verified that on the timescales of interest here, using the GSOD (UTC-based) daily data gives results nearly identical to using the hourly (local-time based) data.
Gaps in the observed temperature record of 5 days or less were filled by linear interpolation, as long as that process resulted in no more than 4 % of the entire data set being filled (in which event the station was rejected). Precipitation had notably more gaps than temperature, and there appeared to be a problem with zero precipitation values occasionally being reported as missing values. Accordingly, for precipitation, gaps of up to 8 days were allowed, and gaps were filled with zero. Stations were included if up to 8 % of the data were infilled. Any station was rejected if, after filling the allowable gaps, more than 10 % of the data remained missing. In the analysis by either water year or season, the individual water year or season was excluded if more than 5 % of that year/season's data was missing after the infilling procedure. Finally, three stations that were located within a few kilometers of other included stations were dropped. Typically, included stations had ∼ 1 % of the temperature data infilled, 2-4 % of the precipitation data infilled, and 0.5-2 % missing data remaining at the end of the infilling procedure. Only 5 of the 74 stations had a significant T dew trend over this time period, close to what would be observed by chance. However, inhomogeneities in the observed record due to instrument changes could affect the long-term trend (Brown and DeGaetano, 2009), an issue not explored here. Figure 3 (left column) shows examples of T dew time series from observations (blue) and estimated by VMS (red) at 4 stations. Ideally, the blue and red curves should coincide. For comparison, T min is also shown (black). With the exception of Tucson (KTUS), which is an arid location, T dew (blue) follows T min (black) quite closely. At Tucson, VMS's implementation of the Kimball-97 correction to T dew clearly improves the simulated values. However at the other stations, VMS's value of T dew (red) is no better than (Billings), or clearly worse than (San Diego and Roberts Field) simply assuming T dew equals T min . This is particularly noticeable at KSAN (San Diego), where the correction severely worsens the estimate of T dew during the summer.
A bispectral analysis of the data (Fig. 3, right column) has some unusual properties, depending on the location. At Billings and Tucson, the squared coherence between observed and modeled T dew is largest at the lowest frequency, which makes sense; the T dew estimation algorithm does better capturing the large seasonal differences than the shortperiod, day to day fluctuations. However at San Diego and Roberts Field, coherence actually drops at the lower frequencies. In other words, the T dew estimation algorithm is doing a worse job at capturing the large amplitude seasonal cycle than it is capturing variability at the 0.1-0.02 day −1 (1/(10-50 day)) time scale. This seems counterintuitive, since a larger swing in T min should provide a more robust base signal to correct with the Kimball-97 algorithm. Finally, it is worth noting that the coherences, although all statistically significant given the approximately 35-yr time series that are available, are rather low values. At San Diego, for instance, no more than 40 % of the variability is captured by the T dew estimation algorithm at any frequency, and at the longer time scales (> 50 days), less than 20 % of the variance is captured. Only in Billings does the explained variance both generally exceed 50 % and reach a peak at the lowest frequencies. This illustrates some of the timescale-dependent deficiencies in the T dew calculation. The mean bias in the estimation of T dew is shown in Fig. 4, and the RMS error in Fig. 5. In summer and autumn the errors are largest along the western coast and smaller in the interior, a geographical distinction also found in the T dew cluster analysis of Robinson (1998) overwhelmingly negative, i.e., VMS tends to estimate a lower dew point (drier atmospheric conditions) than observed, typically by about 2-10 • C, and particularly in the western part of the domain. Although many locations in central and southern California have mean biases ≤ −10 • C, no location has a mean bias of > +5 • C. The largest RMS errors are also concentrated in California, along with a few stations in Oregon and Washington, and are largest in summer and autumn, which include part of the growing season and the bulk of the wildfire season. At San Diego (KSAN) and other Mediterranean climates across California, the poor T dew estimation arises partly due to the way VMS uses a modified version of the Kimball-97 algorithm that implements a 90-day window when computing annualized P . In these locations it is not unusual to have little to no rain for 90 days in the spring and summer, with the result that the VMS-modified algorithm drops the annualized P value to the arbitrary value of 8 cm. Such a low P value results in a far too negative T dew correction during the dry period. The upper left panel of Fig. 3 illustrates this phenomenon. By late June of 1991 (small green arrow on the figure), there had been little precipitation at the station for 90 days, and the corrected T dew (red) discontinuously jumps Hydrol. Earth Syst. Sci., 17, 1833-1850 to a much lower value once the 90-day averaging window advances to cover only the dry period. Averaged across all 74 stations examined here, in March-April-May (MAM) the mean T dew bias (VMS minus observations) for the original Kimball-97 algorithm is 3 % larger than for the VMS modified algorithm. However in the other three seasons, the mean bias is 20-150 % larger using VMS's modified algorithm than found when using the original Kimball-97 algorithm. The simulated values are worse primarily in the Mediterranean coastal climates of central and southern California, locations where Figs. 4 and 5 showed the T dew estimate was especially biased. The biggest effect is in September-October-November (mean error of −2.88 • C for VMS, vs. −1.13 • C for the original algorithm with actual long-term average P ), at the end of the long summer dry period in coastal southern California. Diurnal variation in RH in that season can be very large on the coast, with nighttime moisture recovery with the onshore breeze.

Simulating humidity in warm vs. cold years
Our primary purpose is to understand how VMS might simulate humidity in the future, as the Earth becomes increasingly warmer due to the accumulation of anthropogenic GHGs in the atmosphere. This is distinct from the problems VMS and the Kimball-97 parameterization exhibited in the previous sections, simulating observed historic humidity in western US and coastal Mediterranean climates. There is no way to evaluate future behavior based on historical data alone, but observed data can be used to evaluate whether the VMS humidity algorithm generates systematically different biases depending on the temperature. Ideally, even if the algorithm is biased (Fig. 4), it should have the same bias in warm and cold years. Otherwise, a systematic change in temperatures over time might generate a systematic change in humidity bias over time, which would mimic a humidity trend. This trend could either add to or offset a real-world humidity trend, depending on its sign. Figure 6 shows the mean bias in T dew in the five warmest years at each station, minus the T dew bias in the five coldest years at that station. Water years (October to September) were used in this calculation rather than calendar years since we are most interested in hydrological applications. The entire period used is water years 1976-2009 (34 yr), so this is approximately the difference between the warmest 15 % and coldest 15 % of years. Results are given for the entire water year (left panel), or just the cold (middle panel) and warm (right panel) seasons; the warmest and coldest years are calculated separately for each season. The number of stations with positive and negative differences in bias between the warm and cold years is noted in the lower left corner of each panel in Fig. 6. The results show that during the cold season, the T dew bias is less negative in warm years (p < 0.01). Although the warm season result is in the same sense, it is not significant, and the year-round result shows little difference in the number of stations.
While less bias is desirable, the drawback is that as temperatures increase over time this changing bias could mimic a positive T dew trend. In the cold season, the difference in T dew bias between the five warmest and five coldest years is on the order of 1 • C (Fig. 6). The air temperature difference between the warm and cold years is about 3 • C. So if the western US warms ∼ 3 • C in the cold season as a result of climate change, this may produce a trend of ∼ 1 • C in the simulated T dew . However in the warm season, and especially in the yearly average, this does not seem to be an issue.

CMIP5 global climate model humidity changes
We calculated the annual and seasonal (DJF, MAM, JJA, SON) surface relative humidity (RH) trends in 13 global climate models from the Coupled Model Intercomparison Project, version 5 (CMIP5) data archive (Taylor et al., 2012). The models are listed in Table 1, and include all those with the necessary surface daily humidity, temperature, and precipitation data at the time this work was undertaken. We used   RCP 8.5 (van Vuuren et al., 2011), roughly corresponding to medium and high anthropogenic greenhouse gas emissions scenarios, respectively. The RCP numbers indicate the approximate anthropogenic forcing each scenario experiences, in Watts m −2 , at the end of the century. The period included is 2010-2099. Because of the large number of models and seasons, in this work we show multi-model ensemble averages (MMEAs), which in many measures produces superior results to an individual model, even in a regional context ). Figure 7 shows the seasonal and annual RH trends for the MMEA, along with the number of models that have a negative trend (decreasing RH) in each grid cell. On an annual basis, 9-12 of the 13 models have predominantly negative trends over most of the western US, particularly in the northern part of the domain. RH declines are generally stronger away from the coast. The average trend in spring and sum-mer reaches −0.2 percentage points per decade in RCP 4.5, and −0.8 points per decade in RCP 8.5, with typically 80 % or more of the models agreeing on the sign of the trend. All the models show increasing RH offshore of California in the extreme southwest part of the domain, though the values are small. Several models show increasing summer RH over western Arizona and southern California in RCP 4.5, apparently a monsoonal response, although this trend is not statistically significant in the ensemble average and is reversed in the more strongly forced RCP 8.5 scenario.

VMS-simulated humidity changes
To evaluate VMS's simulation of RH trends, the daily temperature and precipitation data from the 13 global models was bilinearly interpolated to a common 2 × 2 • grid over the western US, applied to VMS, and the difference between the VMS-estimated RH trend and global model RH trend was calculated. VMS's algorithms calculate dew point temperature, which is then converted to vapor pressure. To calculate RH, the variable of interest here, requires an air temperature. It is assumed that RH is calculated with respect to the daily average temperature, T avg . Further, T avg was taken as (T min + T max )/2. This is not exactly true for VMS, which internally calculates temperatures at each hour and then averages them to produce T avg , but is nonetheless a close approximation.
Results are shown in Fig. 8. It is clear that VMS tends to systematically bias the RH trend towards more positive values. Errors are largest (VMS simulating an overly positive RH trend) in the interior western US, roughly over Colorado, Idaho, Montana, and Wyoming, and are about twice as large for the RCP 8.5 forcing scenario as for the RCP 4.5 scenario, which is consistent with the anthropogenic forcing being considerably larger in the latter. The sign of the error is reversed near the coast in many of the models, although coastal errors are generally much smaller than errors in the interior. Errors are weakest in the autumn, and generally stronger and approximately equal in the other three seasons. Along the coast the location of greatest error moves north as the year progresses, peaking in spring in Southern California, summer in Central/Northern California, and autumn in the Pacific Northwest.
These errors can be made more quantitative by computing the histogram of differences between the RH trends at each grid cell, pooled over all the models (Fig. 9). The lefthand column shows the VMS-estimated RH trend minus the original global model trend, i.e., the change in RH trend due to VMS. In addition to imposing notable spread, VMS's humidity algorithms systematically bias the simulated RH trends, with mean biases across the domain and different models of 0.12 percentage points per decade for RCP 4.5 and 0.32 percentage points per decade for RCP 8.5. The stronger the applied GHG forcing, the bigger the positive RH trend bias VMS tends to develop. Considering that the ensemble median RH trend from the full set of CMIP5 global models is on the order of −0.2 to −0.8 percentage points per decade (Fig. 7), these biases are easily large enough to be significant to the simulation's results.
In regional modeling it is common to bias correct and downscale global model data, so the middle column of Fig. 9 shows the effect that bias correcting of the original global model data has on the simulated VMS RH trend (the bias correction method used is described in Pierce et al., 2012). The right column shows the change in trend due to bias correction with constructed analogue (BCCA) downscaling Maurer et al., 2010). For constructing this last quantity, VMS RH trends were simply interpolated to the finer downscaled grid before the grid cell-by-grid cell difference was computed. Neither bias correction nor downscaling contributes appreciably to the error in the RH trend.

Components of the relative humidity trend
Decreasing RH can arise from increasing temperatures or decreasing atmospheric water content. Both T avg and T dew increase in the model runs analyzed here, but T avg always warms more than T dew , typically by about 50 %. So the decrease in RH is accomplished by ambient air temperatures Fig. 9. Histograms of changes (shifts) in the estimated RH trend (percentage points per decade) at all gridpoints in the western US, accumulated across all models. Left column: RH trend in VMS minus that found in the global model. Center column: trend after bias correction minus trend before bias correction. Right column: trend after downscaling minus that found before downscaling. Y-axis is number of grid cells. The mean (x) and standard deviation (σ ) of the distribution are also indicated. warming faster than the increasing moisture content of the atmosphere, but both increase.
Equal trends in T dew and T avg do not contribute equally to a trend in RH because of nonlinearities in the relationship between temperature and humidity. RH = svp(T dew )/svp(T avg ), where svp(T ) is the saturation vapor pressure at temperature T ; since both T dew and T avg have trends, the relative humidity trend dRH/dt can be expanded using the standard formula for the derivative of a quotient, along with the chain rule, to relate RH trends to trends in T dew and T avg . Actually evaluating the terms requires calculating dsvp/dT (note: change with respect to temperature, not time), which was obtained by taking the tangent to the svp(T ) curve at the mean climatological T dew and T avg values (i.e., linearizing around the annual mean T dew and T avg ). Figure 10 shows the resulting contributions to the global model RH trend that arise from global model trends in T dew and T avg . In the RCP 4.5 runs, warming T avg by itself would tend to change RH by −0.6 to −1 percentage points per decade in the majority of the interior western US. T dew by itself would tend to increase RH by +0.5 to +0.8 percentage points per decade, partially compensating for the T avg warming, leaving the residual RH trends

Understanding VMS's humidity trend bias
As described in Sect. 2, VMS computes T dew using T min and a correction based on the local aridity. The first question is whether the humidity error arises from the handling of T min or from the aridity correction. What would the simulated RH trend be if only T min were included in the calculation? We computed this and found that in about half the simulations, using only T min results in an even more positively biased trend than found when including the correction. In the other half, the trend is about the same as found using the correction. The implication is that the correction, while often acting in the right direction, is not strong enough to overcome the deleterious effects of assuming T dew can be estimated simply as T min . This point is brought out more fully in Fig. 11. Panels a-c show the mean global model T min trend, T dew trend, and difference between them, respectively, for the RCP 4.5 forcing scenario. For comparison, the difference between the T min and T dew trends computed by VMS is shown in panel d. In the interior, the difference in trends found in VMS is no-tably smaller than in the global models. In other words, when computing T dew , VMS does not correct the T min trend enough for aridity, and as a result T dew in VMS increases too much in the interior western US. The implied error in VMS's T dew trend is shown in panel e. Similar, though stronger, results are found for RCP 8.5 forcing (panels f-j).
The Kimball-97 aridity correction depends on DTR and EF (Eq. 1). The model-estimated DTR trend is shown in Fig. 12. There is little agreement in the magnitude or sign of this trend across the different models, leading to a nearzero annual mean value for RCP 4.5 and values of 0.05 or less for RCP 8.5. Furthermore, Eq. (1) indicates that a DTR trend of 0.05 • C per decade, consistent with the individual seasonal values of RCP 4.5 or the annual mean for RCP 8.5 as shown in Fig. 12, will produce a T dew tendency of only about 0.01 • C per decade. This value is too small to explain the 0.05-0.1 • C per decade error in the VMS-estimated T dew trend (Fig. 11). Between the model disagreements on the sign of the DTR trend and the small effect it has, it is unlikely that the overly large T dew trends in VMS arise from treatment of the DTR parameter.
The EF trends are also shown in Fig. 12 trends (Fig. 12) to the errors in the VMS-estimated T dew trend (Fig. 11, panels e and j) reveals correspondences; where the EF trend is greatest, the T dew error tends to be lowest. For example, in RCP 4.5, the EF trend is greatest along the west cost and southern tip of Arizona/New Mexico; this is also where the T dew trend error is smallest. The point-by-point spatial correlation between the two fields is 0.66 with a slope that is significantly different from zero (p < 0.05). I.e., locations where VMS simulates increasing aridity (and therefore an increasing depression of T dew below T min ) also show T dew trends that are the most consistent with the original global model. Locations in the interior western US where VMS fails to show increasing aridity show the largest errors in the T dew trend.
This analysis implies that the problem is that not enough of the domain shows the increasing EF (aridity) trends, particularly in the northern interior western US -exactly where VMS's T dew trend is too large. Since EF is PET/P , the low values of P in the southwest deserts drive very large values of EF there. The percentage change in EF is much more equal across the domain than seen in Fig. 12, but the actual magnitudes of the EF changes are much larger in the desert regions because of the division by P . In the Northern Rockies, a different weighting of EF or the EF trend would have little effect on the T dew trend, since the EF trend is so small there to begin with. The problem of insufficient EF trend could be solved by adding a bias term to the EF trend, but that would simply be inserting the anticipated trend into the result in an ad hoc way, and so should be avoided.

Implications for runoff over the western US
The errors in simulated humidity might have an effect on runoff, with lower humidities associated with more evapotranspiration water loss from the surface and therefore lower runoff. To test this, a simple numerical experiment was performed by decreasing VMS's calculated value of T dew by 0.75 • C before it was used by VIC. This value was chosen as a representative end-of-century value for the interior western US based on the mean model analysis (Fig. 11, panels e and j).
For a control run, VMS was first driven by observed daily T min , T max , and P on a 1/8 × 1/8 • latitude-longitude grid over the western US, over the period 1915-2003(Hamlet and Lettenmaier, 2005. The experiment consisted of forcing VMS to decrease its calculated T dew by 0.75 • C at all times, with the intention of isolating the effects of T dew by altering only T dew , leaving T min , T max , and P unchanged.
The effect of decreasing T dew on runoff is shown in Fig. 13. Declines exceed 4 % over much of the western US, ranging up 14 % in some locations. The decreases are not geographically uniform, with the largest values concentrated in the higher elevations east of both the Sierra Nevada and Cascade mountain ranges, and in the mountainous parts of Colorado, Utah, Idaho, Montana, and Arizona. In the lower elevation parts of the lower Colorado River basin, by contrast, decreases generally do not exceed 2 %.
After routing the runoff, annual streamflow in the Colorado River at Lees Ferry declines about 4 %. A recent study (Harding et al., 2012) using VMS forced by downscaled GCM temperature and precipitation fields as described herein concluded that Colorado River streamflow out of the upper basin might decline 7.6 % by the end of the century (a mean value across models), so a 4 percentage point correction to that result is potentially important to include.
Although simplistic, this experiment suggests the humidity errors could have a discernable impact on modelsimulated future streamflow. We emphasize that the concern is not with a constant bias in the simulated Lees Ferry runoff, which could be addressed by a simple bias correction technique. Rather, the issue is that the global models project a trend of decreasing relative humidity in this area in the coming century, a trend that is not well captured by VMS. Correctly including this humidity trend could decrease simulated Colorado River flow by about 4 % over the course of the century.

Summary and conclusions
In this work we have examined relative humidity (RH) simulations from 13 global climate models from the CMIP5 archive, driven with the RCP 4.5 and 8.5 scenarios of atmospheric greenhouse gas concentrations over the course of the 21st century. Results suggest that the interior western US is likely to experience an RH decline of about 0.1 to 0.6 percentage points per decade, depending on season, location, strength of the anthropogenic forcing, and the model considered. Land surface models should include this drying, as it could affect processes such as evapotranspiration, runoff, wildfires, and plant growth. The purpose of this work has been to determine if the variable infiltration capacity (VIC) hydrological modeling system (VMS), when run with input variables T min , T max , and precipitation (P ) from a global model, preserves the original global model's RH trend. We examine this hydrological model because it has been frequently used in studies of the future climate of the western US. Systematic errors in its simulations could therefore have implications for our understanding of climate change's effects in the region, and potentially on resource management or adaptation decisions made on the basis of those studies. VMS calculates atmospheric moisture content using algorithms from the MT-CLIM package (Hungerford et al., 1989;Thornton and Running, 1999;Thornton et al., 2000), incorporating a modified version of the Kimball et al. (1997)   parameterization of T dew . Local T dew is taken as T min , adjusted by an aridity term (arid locations have a lower T dew relative to T min ) and a term proportional to the diurnal temperature range (DTR). Comparison of VMS's estimated T dew to observed values at 74 meteorological stations across the western US shows that the parameterization does a poor job at locations influenced by the marine environment, particularly along the western coast of the US, where it consistently underestimates T dew . This is partly due to neglect of marine humidity sources (and thus not a surprising result), and partly due to a modification to the  algorithm that is used in VMS; rather than the aridity term being calculated from annual precipitation, it is calculated from annualized precipitation over the last 90 days. In the Mediterranean climates along the west coast, it is not unusual to experience 90 days without any precipitation, leading to a sudden drop in VMS-estimated RH. Although these Mediterranean locations typically produce limited runoff (for instance, Southern California imports most of its water), correctly estimating humidity in the population centers there is still of interest for purposes of modeling human comfort (heat index), health, air conditioning use, and wildfire. Wildfire probabilities in neighboring vegetated areas could be misestimated, along with the air pollution impacts of changes in wildfire (both local and regional, since wildfire emissions are transported long distances).
To examine whether a global model's RH trend is preserved in VMS, the daily T min , T max , and P fields from 13 global model simulations were applied to the VIC modeling system. It was found that the VMS-simulated RH trend is not a faithful reproduction of the trend in the original global model, but rather is consistently biased towards positive values in interior locations (i.e., less RH decrease than the global models projected). Since regional climate change studies often use bias correction and downscaling, we also examined the impact these have on the RH trend bias, but found they had negligible effect. The size of the bias is significant, being enough in some cases to eliminate or even reverse the original RH trend found in the global model.
A simple numerical experiment suggested that eliminating the RH bias could diminish runoff up to 14 % in the dry interior US, and reduce Colorado River flow at Lees Ferry by up to 4 % by the end of this century. As of this writing, the Southern Nevada Water Authority (http://www.snwa.com) is constructing infrastructure that will enhance their access to an amount of water equal to 2 % of the historical Colorado River flow at Lees Ferry. That project has an estimated cost of approximately 800 million USD, so we believe that a systematic model bias of 4 % in projections of future Lees Ferry flow is of practical economic significance.
The RH trend errors could also affect the simulation of future wildfires. The fire model in Westerling et al. (2011) finds that the distribution of fire sizes is sensitive to moisture deficit, which is affected by relative humidity. This distribution is highly non-gaussian, with a heavy tail of large events that have a tremendous impact on the western landscape, so a possible change in fire size could be significant for understanding future climate impacts in the western US. In addition to fire size, Westerling et al. (2011) also model large (> 200 ha) fire occurrence as a nonlinear function of moisture deficit. Analyzing their model over the estimation period used in that work , we find that a systematic 5 % decrease in relative humidity from within the range of relative humidity values where fire is more likely to occur increases the average probability of a large fire occurring by 13 to 60 %, with the larger percentage increases corresponding to higher initial relative humidity values (i.e. lower initial probabilities of fire).
Simulated RH declines are slightly too strong along the coast. This is likely the result of the VMS-estimated aridity trends and the corresponding corrections to the T dew − T min relationship failing to represent ocean humidity sources. However, the geographical extent of marine influences may be exaggerated in the relatively coarse resolution global climate models used here, which do not resolve coastal topography.
The global model results agree that both average temperatures (T avg ) and T dew will increase everywhere in the western US in the coming century. Generally T avg increases more rapidly than T dew , leading to a decrease in RH, especially in the interior. VMS estimates an overly large increase in T dew , which is why the RH decrease is smaller in magnitude than it should be.
These results suggest that improving the simulation of RH trends under conditions of climate change may require introducing a new input variable to augment the current dependence on T min , T max , and P . The accumulation of greenhouse gases in the atmosphere is altering the radiation balance of the planet, initially through an increased downward longwave flux from the atmosphere to the surface. A statistical fit to historical conditions might have limitations if the balance of physical processes determining the relationship between T dew and T min changes, and this may be the case for future humidity over the western US.
Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table 1 of this paper) for producing and making available their model output. For CMIP the US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. The authors would also like to thank anonymous reviewers for helpful comments on this manuscript.
Edited by: N. Romano