A new hybrid statistical–dynamical downscaling technique is described to project mid- and end-of-twenty-first-century local precipitation changes associated with 36 global climate models (GCMs) in phase 5 of the Coupled Model Intercomparison Project archive over the greater Los Angeles region. Land-averaged precipitation changes, ensemble-mean changes, and the spread of those changes for both time slices are presented. It is demonstrated that the results are similar to what would be produced if expensive dynamical downscaling techniques were instead applied to all GCMs. Changes in land-averaged ensemble-mean precipitation are near zero for both time slices, reflecting the region’s typical position in the models at the node of oppositely signed large-scale precipitation changes. For both time slices, the intermodel spread of changes is only about 0.2–0.4 times as large as natural interannual variability in the baseline period. A caveat to these conclusions is that interannual variability in the tropical Pacific is generally regarded as a weakness of the GCMs. As a result, there is some chance the GCM responses in the tropical Pacific to a changing climate and associated impacts on Southern California precipitation are not credible. It is subjectively judged that this GCM weakness increases the uncertainty of regional precipitation change, perhaps by as much as 25%. Thus, it cannot be excluded that the possibility that significant regional adaptation challenges related to either a precipitation increase or decrease would arise. However, the most likely downscaled outcome is a small change in local mean precipitation compared to natural variability, with large uncertainty on the sign of the change.
Freshwater in the Los Angeles region comes from local storms, snowpack drainage, and groundwater. Identifying how climate change may impact these sources is of pressing concern for ecosystems and municipal, agricultural, and recreational purposes. In this study, we only aim to quantify twenty-first-century climate change impacts to mean local sources of precipitation across the greater Los Angeles region. Local sources contribute approximately 10% to the water supply in the City of Los Angeles (Villaraigosa 2008). However, in some areas, such as the San Fernando Valley, it contributes a larger portion (Sheng and Wilson 2008; Slade 2013). Furthermore, these local sources may come under increasing pressure in the future (Erb et al. 2011). We do not address potential changes to imported water sources (e.g., the Colorado River) or extreme events (e.g., Das et al. 2013) in this study. A separate study will examine responses of local snowpack to climate change.
Projecting future precipitation changes over the Los Angeles region is challenging for two reasons. First, in GCM projections the region typically lies at the boundary of two oppositely signed, large-scale zones of predicted precipitation change (van Oldenborgh et al. 2014), as described by the “rich get richer” or “wet regions get wetter and dry regions drier” effect (Chou and Neelin 2004; Held and Soden 2006; Trenberth 2011; Durack et al. 2012). Northern, midlatitude areas are projected to get wetter, while southern, subtropical areas are projected to become drier. Second, the complex topography of Southern California creates variations in precipitation that cannot be represented by coarse-resolution GCM simulations. It is particularly important to adequately represent the coastal mountains over Southern California because they generally lead to significant orographic precipitation effects (Hughes et al. 2009; Neiman et al. 2002).
To address the limitations of coarse-resolution GCMs, a common practice is to downscale global projections to a much finer resolution. Dynamical and statistical downscaling techniques are available to perform such a task. Dynamical downscaling solves the equations of motion and other atmospheric equations numerically, using a regional model that is forced along the boundaries by GCM output. This may represent the most physically consistent method to downscale climate data, although systematic biases may still be present in the downscaled simulation. The major trade-off for dynamical downscaling is the expense of huge computational costs. Dynamical downscaling of climate change signals has been done for Southern California. For example, Duffy et al. (2006) dynamically downscaled two GCM projections, finding no statistically significant change in precipitation over Southern California.
Statistical downscaling is computationally cheap compared with dynamical downscaling but hinges on currently existing relationships that may or may not hold true in the future. This technique has also been applied in the region of interest. For example, Hayhoe et al. (2004) statistically downscaled four GCMs using historically derived empirical relationships and found small decreases in future wintertime precipitation in Southern California for three of the four simulations. A recent study by Pierce et al. (2012) uses separate dynamical and statistical downscaling techniques across 16 global climate models to examine future precipitation changes over California. Like Hayhoe et al. (2004), the statistical downscaling approaches used in Pierce et al. (2012) rely only on historical relationships (i.e., they assume stationarity) between variables when calculating climate change signals. After averaging across all downscaled projections, the authors find wintertime precipitation decreases of 5% over Southern California. Maurer (2007) statistically downscale future global precipitation and temperature output to drive a hydrologic model and found slight increases in wintertime precipitation over a basin in Southern California. Das et al. (2013) statistically downscaled 16 GCMs over the Sierra Nevada and found increased 3-day flood discharges, even though models tended to disagree on the sign of mean annual precipitation change. Pierce et al. (2013) also examined possible changes to daily precipitation over California. Using both dynamical and statistical downscaling techniques, they found evidence of increased wintertime precipitation over California, particularly over the northern part, because of an increase in daily precipitation intensity. Note that these previous studies relied on models from phase 3 of the Coupled Model Intercomparison Project (CMIP3), while this study only analyzes models from phase 5 of CMIP (CMIP5). The two ensembles may exhibit different behavior in some cases. For example, Neelin et al. (2013) found that ensemble-mean drying in the CMIP3 archive was stronger over Southern California than in the CMIP5 archive.
The present study uses a new blended dynamical–statistical approach to project mid- and end-of-twenty-first-century December–February (DJFM) precipitation changes at a high resolution over the Los Angeles region. Whereas previous studies use only a dynamical or empirical statistical downscaling technique, this study develops statistical relationships directly from dynamically downscaled output. Using this method we are able to limit the assumption of stationarity that is often employed in statistical downscaling exercises (e.g., Hayhoe et al. 2004; Maurer 2007; Pierce et al. 2012). This technique also allows for downscaling of 36 GCMs in the CMIP5 archive (Table 1), providing analyses on intermodel spread and ensemble-mean changes. In addition to projecting twenty-first-century precipitation changes over Southern California, another major aim of this study is to place climate change signals in context of the region’s significant hydroclimate variability. Huge interannual variability in precipitation over Southern California is largely attributed to its relationships with large-scale natural climate variability patterns such as the El Niño–Southern Oscillation and the Pacific–North American pattern (Cayan and Roads 1984; Redmond and Koch 1991; Dettinger et al. 1998; Cayan et al. 1999; Leung et al. 2003; Berg et al. 2013).
The structure of the study is as follows: Section 2 describes the downscaling techniques and provides observational evaluation of the current climate simulation. Section 3 shows future precipitation changes according to 36 downscaled GCMs and explains the physical mechanisms behind the changes. A discussion of the relationship between climate change and interannual variability patterns is presented in section 4, with a summary of major findings in section 5.
2. Downscaling techniques and validation results
a. Dynamical downscaling
1) Dynamical downscaling framework
A dynamical downscaling simulation over Southern California was performed using the Weather Research and Forecasting Model (WRF), version 3.2 (Skamarock et al. 2008). We use three nested domains (18, 6, and 2 km) to reach a resolution high enough to represent the complex topography and coastlines of Southern California adequately. The three domains and the topography associated with the outermost 18-km domain are presented in Fig. 1a. The outermost domain encompasses all of California and the adjacent Pacific Ocean, while the middle domain focuses on Southern California, including the southern Sierra Nevada. Finally, the innermost, 2-km domain is centered over the greater Los Angeles region. Topography associated with this domain is seen in Fig. 1b.
Model sensitivity experiments were performed to find an optimal WRF configuration. Specifically, this simulation uses the Kain–Fritsch (new Eta) cumulus scheme (Kain 2004), Yonsei University boundary layer scheme (Hong et al. 2006), Purdue–Lin microphysics scheme (Lin et al. 1983), Rapid Radiative Transfer Model longwave radiation (Mlawer et al. 1997), Dudhia shortwave radiation schemes (Dudhia 1989), and Noah land surface model (Chen and Dudhia 2001). Each domain has 43 sigma levels in the vertical and vertical resolution is increased below 3 km to better simulate surface and boundary layer processes.
Two time periods are simulated to initially project mid-twenty-first-century precipitation changes. We focus first on a “baseline” period spanning 1981–2000. In this case, WRF is forced along the boundaries of the outermost domain by the North American Regional Reanalysis (NARR). Then we simulate a range of future climates based on model output from five CMIP5 GCMs (CCSM4, CNRM-CM5, GFDL CM3, MIROC-ESM-CHEM, and MPI-ESM-LR; see expanded model names in Table 1), all under the representative concentration pathway 8.5 (RCP8.5) emissions scenario. For each future simulation, baseline (1981–2000) boundary conditions from NARR are perturbed with future monthly climatological changes (2041–60 average minus 1981–2000 average) to atmospheric variables and imposed on WRF. Three-dimensional atmospheric variables that were perturbed include temperature, relative humidity, zonal and meridional winds, and geopotential heights. Surface temperature, relative humidity, winds, and pressure were also perturbed. This technique has been used previously (e.g., Schär et al. 1996; Hara et al. 2008; Knutson et al. 2008; Kawase et al. 2009; Lauer et al. 2010; Rasmussen et al. 2011; Seo and Xie 2011; Gutmann et al. 2012) and estimates future climates as perturbations to the same baseline mean state, corresponding roughly to the present day. A limitation to this technique is that future interannual variability equals that of the baseline period. Implications of this limitation when analyzing downscaled changes are discussed further in sections 3c and 5. For an application of this downscaling method applied to future warming over the Los Angeles region, the reader is referred to Sun et al. (2014, manuscript submitted to J. Climate).
We first perform a 20-yr future simulation (2041–60), downscaling climate change signals in CCSM4. Computational expenses prevent full 20-yr simulations for other models, so we performed a sensitivity test examining how long of a future period we needed to simulate to capture the full 20-yr climate change signal. Figure 2 shows that by only simulating 3 future years (2058–60) we are able to capture the full 20-yr signal to a high degree of accuracy. (Other consecutive 3-yr periods between 2041 and 2060 may also be highly representative of the full 20-yr time span, although computational resources prevented this analysis.) Spatial structures between the two signals are tightly correlated, with only slight discrepancies seen in the coastal zone. Averaged over the land, the 20- and 3-yr signals are −46.7 and −46.6 mm per wet season, respectively. Relying on this knowledge, we next dynamically downscaled the four other GCMs (CNRM-CM5, GFDL CM3, MIROC-ESM-CHEM, and MPI-ESM-LR) for 2058–60. In each simulation, boundary conditions were created by adding the 2041–60 minus 1981–2000 GCM changes to the 1998–2000 NARR values. Therefore, interannual variability over 2058–60 is the same as 1998–2000; however, the perturbations imposed in the future runs represent a climate change signal associated with much longer averaging periods. Statistical downscaling techniques are then developed based on these 2058–60 minus 1998–2000 dynamically downscaled changes (section 2b). The reader is directed to the supplementary material for a discussion on possible biases in the statistical model due to interannual variability differences between the 3- and 20-yr-long changes. Namely, it is found that the dynamical model underestimates the magnitude of average precipitation changes by around 20% because it is based on 3-yr-long signals (2058–60 − 1998–2000) compared to the 20-yr-long signals (2041–60 − 1981–2000). Thus statistical estimates based on these changes may be associated with a reduced spread by a complimentary amount. The implications of this error are examined further in section 5.
2) Model evaluation: Spatial and temporal variability in the baseline
Before presenting the results of the climate change experiments, we compare simulated interannual precipitation variations in the baseline (1981–2000) 2-km WRF output to observations. We use three observational datasets: California Irrigation Management Information System (CIMIS; http://wwwcimis.water.ca.gov/Default.aspx). National Oceanic and Atmospheric Administration (NOAA) Climate Prediction Center (CPC) 0.25° × 0.25° daily U.S. unified precipitation (http://www.esrl.noaa.gov/psd/data/gridded/data.unified.html), and the 0.5° × 0.5° gridded University of Delaware (UDel) precipitation product (http://www.esrl.noaa.gov/psd/data/gridded/data.UDel_AirT_Precip.html). Correlations between these datasets and WRF output may be less than 1.0 for multiple reasons, including WRF inaccuracies, unresolved subgrid-scale topography (i.e., elevation mismatch between the location being sampled and the WRF gridcell average), poor observational data quality, and inaccuracies in the boundary conditions (NARR) forcing WRF. Assuming the observational products are perfect, the model evaluation serves as a test of WRF’s ability to reproduce precipitation variations over the Los Angeles region when coarse-resolution conditions (NARR) are imposed on it. If WRF is able to transform this coarse-resolution data into regional climate information that closely matches accurate observational products, we are confident WRF can regionalize the GCM signal in a way that is consistent with the real atmosphere’s dynamics.
In Fig. 3a, we correlate monthly DJFM precipitation accumulations in the baseline period between each CIMIS station and the nearest WRF grid point. Each correlation in based on a maximum sample size of 80 (4 wet-season months × 20 baseline years = 80 values). However, there are missing values in the observations, leading to an average sample size of 45 values. Of the 13 stations, 12 have correlations to WRF above 0.5 and more than half have correlations above 0.7. Thus, WRF generally simulates monthly precipitation variations at rain gauges across the domain reasonably well. The lone exception is Santa Barbara (r = 0.37). We speculate that WRF simulates the complex interactions between small-scale circulations and rainfall at this location of intense coastal topography poorly. In Fig. 3b, we correlate 1981–2000 DJFM-mean precipitation accumulations (20 values per grid point) between each CPC grid point and the nearest corresponding WRF grid point. Correlations greater than 0.6 are found across nearly the entire domain, with very high values (r > 0.9) found along much of the densely populated coastal region. The domain-average correlation is 0.82. Thus, the interannual variability simulated in WRF and that recorded in the CPC gridded product are very similar.
Additional validation of precipitation variability in the baseline WRF simulation is presented in Fig. 4. This figure compares interannual variability of monthly precipitation amounts in the three observational datasets (CIMIS, CPC, and UDel) and WRF output at the scale of the domain. Each white, gray, or black dot in Fig. 4 represents monthly precipitation accumulations for each of the 20 baseline years that are simulated. The large dots represent monthly climatologies for each dataset. Two comparisons can be made in Fig. 4. The first is between CIMIS station-averaged monthly precipitation accumulations (white dots; see Fig. 3a for station locations) and corresponding accumulations averaged over the nearest grid points in the 2-km WRF domain (light gray dots). The levels of interannual variability in CIMIS and WRF station averages are very similar for each month, and the two time series are highly correlated (r = 0.88). Climatological accumulations for each month are also very similar, with an average monthly climatology difference between the two datasets of approximately 6 mm, or 8%. Particularly noteworthy is the similarity between the observed and modeled bimodal structure of the temporal precipitation distribution, seen most dramatically in January and February. Both datasets capture the extremely dry (<25 mm) and wet (>250 mm) months within the baseline period.
The second comparison to make in Fig. 4 is between the UDel, CPC, and WRF land-average monthly accumulations (medium gray, dark gray, and black dots, respectively). Like the CIMIS comparison, WRF variability in monthly precipitation accumulations tightly matches what is observed in the UDel (average r = 0.94) and CPC (average r = 0.96) datasets. Differences in monthly climatologies between WRF and UDel are approximately 17 mm (28%) and approximately 9 mm (15%) between WRF and CPC. Interestingly, for both WRF-based and observation-based datasets, there are strong similarities in magnitude between the station-averaged (white and light gray dots) and land-averaged values (medium gray, dark gray, and black dots). This indicates that the station-averages adequately sample the land fraction of the domain. For example, the average monthly climatology difference between CIMIS station-averaged (white dots) and CPC land-averaged (dark gray dots) values is only approximately 16 mm.
Finally, we assess WRF’s ability to simulate spatial variations in station-averaged (in the case of CIMIS rain gauges) and land-averaged (in the case of UDel and CPC gridded observations) precipitation totals over the baseline period. Results are seen in Fig. 5, which shows scatterplots between simulated and observed (CIMIS: black circles; UDel: red circles; and CPC: cyan circles) station or land-averaged wet-season total accumulations. Note that CIMIS observations begin in 1989, so only 12 wet seasons are included in this portion of the plot. WRF reproduces the CIMIS observations (r = 0.83; average bias of +15 mm) better than UDel (r = 0.59; bias of +229 mm) or CPC (r = 0.55; bias of +221 mm). The large disagreement between WRF and the two gridded products is likely due to the horizontal-resolution differences between them. Coarse resolutions in the gridded products (0.25° × 0.25° for CPC and 0.5° × 0.5° for UDel) may not resolve the full orographic effects on precipitation, which are included in WRF and of course the station measurements. As noted above, discrepancies between WRF and CIMIS values or any data product may arise because of subgrid-scale topography and poor observational data quality, in addition to model deficiencies and inaccuracies in the boundary conditions (NARR).
b. Hybrid dynamical–statistical downscaling framework
1) Empirical orthogonal function analysis
Here we present the hybrid dynamical–statistical approach to generating future precipitation projections. We begin by forming statistical relationships between precipitation changes in the five dynamically downscaled GCMs to large-scale parameters in GCM output. The first step is identifying common spatial patterns between monthly wet-season precipitation changes (2058–60 minus 1998–2000) for all five downscaled models. Each GCM’s dynamically downscaled monthly precipitation changes over the course of the wet season (DJFM) can be seen in Fig. 6. We make two remarks on the variations in Fig. 6. First, there is variation in the sign and magnitude of mid-twenty-first-century precipitation changes in dynamically downscaled results. Some downscaled GCMs, such as CCSM4 (Fig. 6, top), show future drying over most of the coastal zone and high elevations for all months, while others, such as CNRM-CM5 (Fig. 6, second row from top), project moistening for much of the domain over most months. Other outcomes lie between these two cases, and are not necessarily consistent in sign across the domain. Second, we note that, although there is large variation across downscaled models and months, there appears to be a common area where most of the action occurs: a pattern tied to orography, with enhanced loading in the coastal zone and throughout the mountainous regions. This suggests that performing an empirical orthogonal function (EOF) analysis on the aggregated set of these monthly precipitation change patterns could yield a single, robust spatial pattern of change.
Following this reasoning, an EOF analysis is performed over the spatial patterns in Fig. 6. Since the EOF analysis spans both models and months, the patterns it generates maximize both intermodel and intermonthly variability. The three leading modes are shown in Fig. 7. The first accounts for 70% of the variability seen in Fig. 6, confirming our suspicion that the majority of the variance can be accounted for with a single spatial pattern. Mode 1 physically represents the dominant orographic pattern of precipitation over Southern California (Hughes et al. 2009; Conil and Hall 2006; Neiman et al. 2002). Other precipitation patterns, such as “blocked events” (Hughes et al. 2009; Neiman et al. 2002), may be somewhat represented in modes 2 and 3. A corresponding 20-value (5 dynamically downscaled models × 4 months) series of mode 1 loadings is also produced from the EOF analysis. These loadings represent the contribution of the spatial pattern of mode 1 to each model’s monthly precipitation change. Since this mode accounts for the majority of intermodel and intermonthly variability, it should be possible to predict the dynamically downscaled precipitation changes in Fig. 6 with reasonable accuracy simply by multiplying the spatial pattern of mode 1 by each model’s monthly mode 1 loading. (While modes 2 and 3 may represent a physical phenomenon associated with precipitation change, we ignore them because of the small variance that is captured in each mode, 7% and 5%, respectively.) Blending the statistical methods of an EOF analysis and dynamical downscaled simulations forms what we call a hybrid dynamical–statistical downscaling technique. For an example of how this blended statistical–dynamical downscaling approach can be applied to regional warming patterns, the reader is referred to Walton et al. (2014, manuscript submitted to J. Climate).
2) Predicting mode 1 loadings
We have calculated mode 1 loadings for the five dynamically downscaled models, but we need a method for predicting the mode 1 loadings for the other GCMs if they were dynamically downscaled. The first step is to relate the known mode 1 loadings to a large-scale predictor variable available from the GCMs, in this case precipitation. In Fig. 8a, we correlate mid-twenty-first-century monthly DJFM precipitation changes over the North Pacific in the five GCMs that were dynamically downscaled to the loading series associated with mode 1. Each GCM is regridded to a common horizontal resolution (1.5° × 1.5°) before performing the correlation. A dipole correlation pattern emerges. GCM precipitation change over the Gulf of Alaska shows anticorrelations to regional precipitation changes associated with mode 1, while the Pacific Ocean adjacent to California shows positive correlations. A physical interpretation of this correlation pattern is discussed in section 4b. We tried several statistical techniques to relate mode 1 loadings to GCM precipitation changes, including single and multivariable linear regression and a projection-based dot product technique. The strongest (highest correlations) and most robust (across the downscaled GCMs) relationship was found using linear regression, where mode 1 loadings are predicted by two predictor variables: GCM precipitation changes averaged over the two regions spanning the dipole correlation pattern (black boxes in Fig. 8a). This yields a single equation to predict a given GCM’s mode 1 loading, if that GCM were dynamically downscaled, based only on its mid-twenty-first-century precipitation change across the northeast Pacific Ocean. A caveat is that these predictive equations hinge on the training set of dynamically downscaled models: in this case CCSM4, CNRM-CM5, GFDL CM3, MIROC-ESM-CHEM, and MPI-ESM-LR. A different set of models could give different relationships between GCM and local precipitation changes. However, some robustness is provided to the predictive relationships by developing them on a set of GCMs that span the dry-to-wet parameter space (see Figs. 6 and 10).
3) Validating statistical downscaling techniques
The statistical model may capture dynamical model output imperfectly for two reasons: 1) mode 1 is an imperfect representation of regional precipitation change, and 2) it is impossible to predict mode 1 loadings perfectly. Knowing the loadings associated with mode 1 from our EOF analysis of dynamically downscaled simulations, we can test how accurate DJFM-mean changes are based solely on mode 1: that is, the first source of error. This comparison is shown in Fig. 9. Recall that the EOF analysis is performed over monthly changes, so DJFM-mean values shown here are calculated by averaging individual monthly patterns to produce a seasonal mean. First we compare the spatial patterns between the dynamically downscaled changes (Fig. 9a) and those based on mode 1 (Fig. 9b). In general the spatial patterns are very well correlated, aside from modest discrepancies in the Mojave Desert regions. WRF (y axis) versus mode 1–based (x axis) precipitation changes from Figs. 9a,b, now averaged over land, are scattered in Fig. 9c. Mode 1 captures the land-averaged precipitation change extremely well, with the mode 1 changes and the WRF changes falling almost perfectly on the line y = x. These results confirm that, if we have perfect knowledge of mode 1 loadings, then statistically downscaled ensemble-mean changes and the spread in these changes are highly representative of the corresponding dynamically downscaled changes. This is especially true when considering the change averaged over the region’s land areas.
Next we analyze the errors associated with imperfect predictions of mode 1 loadings (i.e., the second source of error in the statistical model) using cross-validation experiments. These experiments use differing subsets of the five dynamically downscaled output to develop a predictive equation for mode 1 loadings. We then predict mode 1 loadings for all dynamically downscaled models and compare them to the actual loadings. Specifically, we perform five experiments. The experiment number is equal to the number of dynamically downscaled models used to determine mode 1 loadings. Each experiment is performed for a varying number of trial runs, consistent with the number of ways it is possible to combine the models. For example, experiment 1 uses one model set of DJFM monthly precipitation changes to determine mode 1 loadings (i.e., any one row in Fig. 6). It has five trials since there are five possible DJFM monthly change values that can be used to predict mode 1 loadings. Experiment 2 uses two model sets of DJFM monthly changes to predict mode 1 loadings for all models, yielding 10 unique combinations (i.e., any two rows in Fig. 6). Experiments 3 (i.e., any three rows in Fig. 6) and 4 (i.e., any four rows in Fig. 6) have 10 and 5 trials, respectively, and experiment 5 (all rows in Fig. 6) has only 1 trial.
In essence, we are testing the robustness of the statistical model as more and more dynamically downscaled information is included in its training. For each trial run in each experiment, we perform all analyses described in section 2b(2) for the dynamically downscaled models being used for mode 1 predictions. That is, we first perform an EOF analysis over the spatial patterns of monthly precipitation changes (e.g., 4 patterns per trial in experiment 1). The EOF analysis yields a series of mode 1 loadings, which are then correlated to the corresponding GCM mid-twenty-first-century precipitation changes across the Pacific Ocean. Finally, GCM mid-twenty-first-century precipitation changes over the regions of maximum positive and negative correlation (which varies according to each trial’s correlation map but is similar to Fig. 8a for all trials) are regressed against that trial’s mode 1 loadings. This yields a predictive equation for mode 1 loadings for each of the five dynamically downscaled models, which can be compared to the known mode 1 loadings.
Table 2 summarizes the uncertainty of the statistical model due to errors in the predictions of mode 1 loadings. The error averaged over all models for all trials is shown in the right column. Errors decrease steadily as the number of models used in the EOF analysis increases. This makes sense, since more intermonthly, intermodel variability is included as more information is fed into the analyses. Specifically, average error is reduced from over 100% when using just one or two models, to just 13% when using five models. It can appear that using four models gives a smaller percent error (−2%) than when five models are used (−13%). However, this simply reflects an average over a very large range when four models are used (from −483% to +466%) compared to a much smaller range using five models (from −103% to +137%).
4) Value added over bilinear interpolation
Here we justify the development of our hybrid statistical–dynamical downscaling technique by comparing results to a simple bilinear regression of the raw GCM data down to 2 km. Figure 10 provides evidence that the hybrid downscaling technique adds significant value in spatial patterns compared to bilinearly interpolating GCM data over Southern California. For each GCM in Fig. 10, spatial patterns that emerge in the interpolated results are broad in scale and have no way of capturing the leading spatial pattern seen in the dynamical downscaling associated with orographic effects. Land-averaged changes between the interpolated and dynamically downscaled GCMs can be quite similar and in some cases closer (likely fortuitous given the coarse GCM resolution) than the corresponding statistically downscaled changes. However, orographic influences on precipitation (e.g., Hughes et al. 2009) are simply not captured in either the raw or interpolated GCM data. Conversely, the hybrid dynamical–statistical downscaling technique is able to capture the orographic imprint on precipitation changes with reasonable accuracy. It should also be noted that the standard deviation between the statistically and dynamically downscaled land-averaged changes is 6.5 and 9.5 mm per wet season, respectively. Thus, the statistical model may underestimate the spread of changes on the order of 30%. We will assess the implications of this potential error in section 4a.
3. Statistical–dynamical downscaling results
Here we predict the regional precipitation projections for all 36 GCMs (Table 1), using the statistical model described in the previous section.
a. Mid-twenty-first-century changes
Mid-twenty-first-century DJFM-mean precipitation changes from all 36 downscaled GCMs are shown in Fig. 11. Recall that the downscaled projections in Fig. 11 are forced to have the same spatial pattern (that of mode 1; Fig. 7) and that the spatial pattern is dialed up or down based on the predicted loading for that GCM. Precipitation changes projected using full dynamical downscaling would have somewhat more spatial heterogeneity than those shown in Fig. 11. Thus, we do not focus on the spatial patterns of change, but rather interpret results from a land-average perspective. The land average can be predicted by the statistical model with a high degree of accuracy once mode 1 loadings are known (see section 2b).
Figure 11 shows an apparently large range of projected changes across models: 13 models project increased precipitation (average of +8.1 mm per wet season) and 23 models project decreased precipitation (average of −8.5 mm per wet season). The most extreme models are MIROC5 and IPSL-CM5A-MR, which project changes of approximately +19 and −25 mm per wet season across the land, respectively. The ensemble-mean land-average change is −2.5 mm per wet season, reflecting a large degree of cancellation between moistening and drying tendencies. Note that the statistical model may underestimate the spread of changes up to 50% (see section 2 of supplementary material), so the true intermodel variability of changes may be 50% larger than described here.
b. End-of-twenty-first-century changes
The statistical model can also be used to project end-of-century (2081–2100 − 1981–2000) precipitation changes. As seen by the dark blue dots in Fig. 4, the ensemble-mean change is near zero for each month and the spread of those changes is smaller than current levels of variability, similar to the midcentury case. In addition to downscaled changes, we also present interpolated GCM changes in Fig. 4 (light blue dots). Like the mid-twenty-first-century changes, the ensemble-mean change by the end of the twenty-first century is near zero for each month. Taken as a whole, Fig. 4 indicates that the model-average downscaled and interpolated GCM scenario for the Los Angeles region is very little precipitation change throughout the twenty-first century.
c. Physical mechanisms
Recall that the dynamically downscaled future simulations are based on GCM perturbations to NARR climatologies along the lateral boundaries for three-dimensional temperature, moisture, winds, and geopotential heights and for surface winds, moisture, and pressure. Interpreting the physical mechanisms behind precipitation changes are constrained to this methodology. While storms entering the outermost domain in WRF (Fig. 1a) are structurally identical between the baseline and future simulations [a limitation raised in Rasmussen et al. (2011)], the perturbation method allow for storms to evolve differently as they propagate toward the innermost domain. For example, possible changes to large-scale circulations, such as the jet stream, can be captured in the three-dimensional wind and geopotential height perturbations. Moisture content in future storms could also change because of perturbations in relative humidity. Future storm strength could also be modified by perturbations to surface pressure and three-dimensional geopotential heights. Bearing this information in mind, we revisit Fig. 8 to identify the physical mechanisms underpinning downscaled precipitation changes (Fig. 11).
As described in section 2b(2), Fig. 8a shows that precipitation changes over Los Angeles are related to large-scale precipitation changes over extreme northern and north-central portions of the eastern Pacific Ocean. The patterns in Fig. 8a suggest that average jet stream position changes across the Pacific Ocean are largely controlling precipitation changes over Los Angeles. A recent study by Neelin et al. (2013) analyzed the relationship between end-of-century California December–February (DJF) precipitation changes and 200-hPa zonal wind speed changes over the northeast Pacific Ocean in 15 CMIP5 GCMs (cf. Fig. 1 of Neelin et al. 2013). Precipitation changes over the California land–ocean region are found to be significantly related to changes in the jet stream (i.e., 200-hPa zonal winds) and associated storm tracks. GCMs projecting increased jet stream wind speeds, associated with an eastward and poleward jet extension, tend to steer more storms toward the coast and lead to overall precipitation increases in this region. GCMs that show weak eastward jet extension and/or wind speed enhancement are associated with minimal precipitation changes. Specifically, the authors find a correlation of 0.76 between end-of-century DJF precipitation changes over California and 200-hPa zonal wind speed over a certain region of the northwest Pacific.
Although our domain of interest is the Los Angeles region rather than the whole state of California, we follow the arguments presented in Neelin et al. (2013) and perform an analysis relating GCM 200-hPa zonal wind speed changes to downscaled precipitation changes. The 200-hPa zonal wind speed changes (2041–60 minus 1981–2000) for the 36 downscaled models are correlated at each grid point in the GCM domain to the domain-averaged downscaled precipitation changes. Each GCM is regridded to a common horizontal resolution (1.5° × 1.5°) before performing the correlation. The results are shown in Fig. 8b. Strong negative correlations are seen across most of the Gulf of Alaska and into western Canada. Conversely, strong positive correlations are seen across the entire north-central Pacific Ocean, centered on Hawaii. This dipole pattern echoes the results found in Neelin et al. (2013) and indicates how jet stream positioning and strength influence future precipitation over the Los Angeles region. Specifically, GCMs that project regional increases (decreases) in jet stream strength off the coast of Southern California lead to increased (decreased) precipitation over Los Angeles.
4. Connection to interannual variability
a. Context of current interannual variability
Here we place the intermodel spread of future precipitation changes in the context of the region’s natural precipitation variability. Examining Fig. 4, we compare the variability across statistically downscaled model projections of future changes (red dots) and levels of interannual variability for the wet season (black dots). Averaged across each month, the standard deviations for the downscaled midcentury precipitation changes are 15, 15, 12, and 14 mm per wet season, respectively. (The standard deviations of end-of-century values are very similar.) The standard deviation of baseline interannual variability of WRF land-averaged monthly averaged accumulations (black dots in Fig. 4) is 61 mm per wet season. Thus, the intermodel variations of downscaled future changes in average precipitation are roughly 25% of the current interannual variability. As noted in section 2, the statistical model may underestimate the standard deviation of the precipitation changes, because of imperfect knowledge of mode 1 loadings, probably by about 30%. So potentially the true standard deviation of precipitation changes is roughly 40% of the variability. However, even after factoring in this possible bias, it is clear that the interannual precipitation variability is large compared to potential changes in the mean. Of course, the mean changes would be sustained on time scales much longer than 1 yr, potentially leading to adaptation challenges. For example, the downscaled models with the most extreme drying and moistening tendencies are associated with mean precipitation changes on the order of 10%. However, such challenges would only materialize if the more extreme models are correct; the average downscaled and interpolated GCM outcome is virtually no precipitation change for the entire century.
b. Relationship between future climate changes and interannual variability
So far we have argued that GCM placement of jet stream and storm tracks in the North Pacific Ocean is the main driver of intermodel variability in future precipitation changes over Los Angeles. Previous studies have also shown jet stream placement, strength, and storm-track steering over the Pacific Ocean can shift because of natural climate variability patterns (Chen and van den Dool 1997; Straus and Shukla 1997; Held et al. 1989). These jet stream and storm-track shifts impact the amount of precipitation over Southern California (Berg et al. 2013; Athanasiadis et al. 2010). The importance of the jet stream for future precipitation change suggests a tight link between the physical underpinnings of interannual variability and simulated climate change.
We begin addressing the relationship between interannual and intermodel variability by analyzing baseline DJFM precipitation from the 1981–2000 WRF simulation forced by NARR. An EOF analysis is performed over 20 spatial patterns of DJFM-averaged precipitation anomalies corresponding to each year of the baseline simulation. The patterns are calculated as anomalies relative to the 1981–2000 DJFM climatology. The leading mode accounts for 86% of the variability, and the corresponding spatial pattern is very similar to the first mode of intermodel variability determined from the climate change experiments (Fig. 12). The leading modes of variability in both the baseline and future cases reflect the strong orographic enhancement of precipitation and the influence of blocking in the coastal zone across the greater Los Angeles region (Hughes et al. 2009). After performing the EOF analysis over the baseline precipitation fields, we then correlate the time series associated with mode 1 (Fig. 12a) to 1981–2000 precipitation anomalies at each grid point in the NARR data. These correlation coefficients are plotted in Fig. 8c and can be compared to the future case (Fig. 8a; section 3). Both cases show a tongue of positive correlations that extend from the coast of California westward into the Pacific Ocean. This tongue is then flanked on the north and south by large swaths of anticorrelations. We also perform a correlation between baseline precipitation and 200-hPa zonal wind anomalies in the NARR data (Fig. 8d) and compare it to the corresponding case associated with future changes in the GCMs (Fig. 8b; section 3c). Both cases show a dipole pattern of large positive correlations across the southern half of the eastern Pacific Ocean and large negative correlations in the northern half.
Such similarities in Fig. 8 confirm that the dynamics of baseline interannual variability are nearly identical to those underpinning future intermodel uncertainty. That is, the region’s precipitation currently vacillates between wet and dry periods with a pattern heavily modulated by orography. The vacillations are largely due to natural variations in the position and strength of the jet stream and subsequent storm-track steering. Models that tend to deflect the jet stream and storms away from Southern California yield drier climates in the future, while models showing a tendency toward jet stream strengthening and increased storm activity over Southern California project a wetter climate. Thus the collection of moistening and drying tendencies in the CMIP5 ensemble can likely be understood as an “excitation” of a natural mode of variability.
5. Concluding remarks
This study uses a hybrid dynamical–statistical downscaling technique to examine mid- and end-of-twenty-first-century precipitation changes over the greater Los Angeles region under the RCP8.5 emissions scenario. Modeling dynamically downscaled precipitation changes with statistical methods, we downscale 36 GCMs in the CMIP5 archive based on changes in each model’s large-scale precipitation fields. There are three major findings of this study. First, the ensemble-mean change for both time slices is essentially zero. Second, while downscaled CMIP5 models disagree on both the sign and magnitude of future precipitation changes over Los Angeles, the spread of possible changes is modest compared to current levels of variability. For both time slices, the statistical model estimates that the standard deviation of land-averaged precipitation change is about 0.2–0.25 of the standard deviation of the interannual variability. As shown in section 2, the statistical model may underestimate the intermodel spread by as much as 30% because of imperfect knowledge of mode 1 loadings. So the true standard deviation of the precipitation change, if all GCMs were downscaled dynamically, could be closer to 0.4 of the interannual variability standard deviation. Thus, even after allowing for potential error in the statistical model, current shifts between wet and dry years are greater than average changes in even the most extreme model projections. However, the sustained moistening or drying seen in the most extreme statistically downscaled models could lead to adaptation challenges. Although these changes are unlikely, they amount to roughly 10% changes in mean precipitation for both time slices. Finally, robust similarities are found between the intermodel variability of future changes and interannual variability of baseline precipitation anomalies. Jet stream placement and strength currently dictates winter precipitation amounts, and also dictates the sign and magnitude of future precipitation changes. To the degree there is uncertainty in future precipitation change over the Los Angeles region, it is because of differences in the simulated response of this phenomenon to anthropogenic forcing.
While there is a great opportunity to assign probabilities of future changes based on an ensemble of projection outcomes, no single method perfectly accomplishes this task. Ensemble-mean forecasts have proven skillful in producing most likely outcomes for climate variations on the seasonal time scale (e.g., Doblas-Reyes et al. 2003; Palmer et al. 2005), along with hurricane paths and other weather patterns (e.g., Zhang and Krishnamurti 1999; Krishnamurti et al. 2000). Interpreting an ensemble-mean projection as the most likely outcome for twenty-first-century precipitation, however, is complicated due in part to the fact that the range of outcomes in the ensemble spans positive and negative changes (van Oldenborgh et al. 2014, their Fig. A1.6). Our result of near-zero ensemble-mean change simply reflects offsetting tendencies of moistening and drying in the GCMs. It could be misleading not to acknowledge that some change, either positive or negative, is likely to occur. At the same time, for nearly all projections, the magnitude of the change is small compared to natural variability. In this sense, we interpret the most likely scenario as a small change in precipitation compared to natural variability, with large uncertainty on the sign of the change.
A critique to this probabilistic estimate is that we are weighting each projection equally in the ensemble, assuming the quality of each model is the same (Tebaldi and Knutti 2007; Knutti 2010; Shepherd 2014). One way to address this issue is to weight each model based on its ability to simulate historical precipitation or the large-scale circulations controlling precipitation over California. Swain et al. (2015) identified 12 CMIP5 GCMs that best simulate historical 500-hPa geopotential height over the northeastern Pacific Ocean according to reanalysis distributions. Geopotential height in this region strongly influences California’s precipitation (Fig. 8). By examining downscaled projections from these 12 “good” models in Fig. 11, we again find that precipitation changes are modest compared to natural variability and that the sign of the change is uncertain: specifically, 8 models project increased precipitation, 4 models project decreased precipitation, and the ensemble-mean change is +2.2 mm per wet season. [Similar results between the 12- and 36-model ensembles are perhaps not too surprising as Langford et al. (2014) found that most CMIP5 models have reasonable skill in simulating California precipitation.] This simple analysis indicates that the major conclusions of this study are not sensitive to which models are included in the ensemble.
Our result of near-zero ensemble-mean precipitation change over Los Angeles can be interpreted in terms of the well-accepted understanding of global precipitation change whereby patterns of precipitation become enhanced, such that wet regions become wetter and dry regions become drier (Chou and Neelin 2004; Neelin et al. 2006; Held and Soden 2006). This leads to increased precipitation over convection zones and drying outside of the convection zones. On average, Southern California is positioned between areas dominated by these competing tendencies: increased precipitation to its north in the midlatitudes and decreased precipitation to the south within the subtropics. However, in some GCMs the region is north of the boundary between the two zones, while in others it is south of it. As such, precipitation projections over this region tend to negate one another and yield small ensemble-mean projections.
One interesting finding from this study is that intermodel variability between the statistically downscaled (red dots in Fig. 4) changes is approximately half the size of the variability according to the GCM-interpolated changes (pink dots). We attribute this spread suppression in the statistical model to two sources. First, inaccurate mode 1 predictions by the statistical model (section 2) can underestimate the intermodel spread by at least 30%. Second, statistical relationships are derived from 3-yr-long dynamically downscaled changes (2058–60 − 1998–2000). On average, 1998–2000 was found to have somewhat lower precipitation than 1981–2000. This has the effect of reducing the magnitude of statistically estimated climate change signals, and the spread associated with them, likely about 20% (see Fig. S1 and Table S1 of the supplementary material). If these two error sources were eliminated, the spread would probably increase by about 50%: therefore, approximately accounting for the difference between statistically downscaled and GCM-interpolated intermodel variability seen in Fig. 4. We also note that the results in this study are limited to our choice of WRF as the regional model. Unlike larger regional modeling efforts such as the Coordinated Regional Climate Downscaling Experiment (CORDEX; e.g., Nikulin et al. 2012) and North American Regional Climate Change Assessment Program (NARCCAP; e.g., Wang et al. 2009), computational costs prevented an examination of how sensitive this study’s results are to multiple regional model simulations. However, we found strong agreement between our results and those from an independent study that statistically downscaled precipitation changes over the greater Los Angeles region (L. Alexanderson, Los Angeles County Department of Public Works, 2014, personal communication).
Differences between the regional model outcomes and those of the GCMs may also stem from our method of perturbing baseline boundary conditions using future climatological changes. For example, one could instead directly downscale raw historical and future GCM data to calculate changes, as opposed to perturbing baseline conditions derived from reanalysis. We are currently conducting research to test whether this direct method gives different results from downscaling changes in the climatology through a perturbation to reanalysis-based boundary conditions.
Given the agreement between the GCMs and the downscaled information in the ensemble-mean outcome, it seems unlikely that a different dynamical downscaling technique would generate a systemically different answer. The hybrid statistical–dynamical downscaling technique could be applied beyond the Los Angeles region. It may be especially appropriate in areas that share these two characteristics with the domain of interest in our study: 1) changes in the large-scale circulation govern precipitation change, allowing for development of credible GCM scaling factors, and 2) local precipitation changes are heavily influenced by orography, leading to diagnosed local response patterns, as encapsulated by the leading EOF patterns. Thus, it would be applicable for any mid-to-high-latitude location with significant topography.
An important caveat relating to the El Niño–Southern Oscillation (ENSO) phenomenon applies to the conclusions of this study. In the current climate, ENSO influences the position of the Northern Hemisphere jet stream and storm tracks across the eastern Pacific Ocean through atmospheric teleconnections (Held et al. 1989; Chen and van den Dool 1997; Straus and Shukla 1997). These shifts have a statistically detectable effect on precipitation over Southern California. During La Niña events, the jet tends to move northward toward the Gulf of Alaska, leading to drier than average conditions across Southern California. Under El Niño conditions, the jet tends to extend south and eastward, steering storms more directly across southern regions of the United States, including Southern California (Redmond and Koch 1991; Dettinger et al. 1998; Cayan et al. 1999; Leung et al. 2003; Berg et al. 2013). The CMIP5 ensemble of GCMs has shown improvements in the simulation of ENSO compared to the CMIP3 ensemble, particularly in the amplitude and time scale of the phenomenon. However, the CMIP5 models still exhibit significant errors, especially in the irregularity of the phenomenon and its spatial pattern (Flato et al. 2014). A detailed examination of the implications of these tropical Pacific errors for precipitation change over Southern California is beyond the scope of this study. However, it seems possible that the GCM projections of future ENSO behavior may be affected by them.
If these errors were corrected, modestly different outcomes for Southern California precipitation might result, because of the link between ENSO variability and Southern California precipitation. When an ENSO event occurs, it accounts for roughly two-thirds of the variance in Southern California precipitation. However, only about 40% of wet seasons can be considered strong ENSO events (Schonher and Nicholson 1989). Thus, roughly one-quarter of the variance of Southern California precipitation can be traced to ENSO. The remaining three-quarters of the variance is linked to shifts of the jet stream unrelated to tropical Pacific variability, similar to those portrayed in Fig. 8d, which are also the mechanism generating intermodal spread in the CMIP5 ensemble. While ENSO is a mechanism generating regional precipitation variability, it is not the most important. ENSO errors in the GCMs may introduce somewhat more uncertainty in our regional precipitation projections than what is implied by the downscaled intermodel spread alone. It is impossible to quantify this effect precisely with present knowledge, but the role ENSO currently plays in Southern California precipitation does at least offer a useful guide. We estimate that ENSO GCM errors increase the uncertainty by an amount approximately proportional to the fraction of the variance ENSO accounts for in current climate: by about 25%. This additional uncertainty underscores the need for regional planning that allows for a variety of future precipitation change outcomes.
Support for this work was provided by the City of Los Angeles and the U.S. Department of Energy as part of the American Recovery and Reinvestment Act of 2009. Additional funding was provided by the National Science Foundation (Grant EF-1065863, “Collaborative research: Do microenvironments govern macroecology?”) and the Southwest Climate Science Center. Additional support in part by National Science Foundation Grant AGS-1102838. We thank two anonymous reviewers for their constructive comments.