A Method for Detection and Attribution of Regional Precipitation Change Using Granger Causality Application to the United States Historical Record

Despite the emerging influence of anthropogenic climate change on the global water cycle, at regional 1 scales the combination of observational uncertainty, large internal variability, and modeling uncertainty undermine 2 robust statements regarding the human influence on precipitation. Here, we propose a novel approach to regional 3 detection and attribution (D&A) for precipitation, starting with the contiguous United States (CONUS) where ob4 servational uncertainty is minimized. In a single framework, we simultaneously detect systematic trends in mean 5 and extreme precipitation, attribute trends to anthropogenic forcings, compute the effects of forcings as a function 6 of time, and map the effects of individual forcings. We use output from global climate models in a perfect-data 7 sense to conduct a set of tests that yield a parsimonious representation for characterizing seasonal precipitation 8 over the CONUS for the historical record (1900 to present day). In doing so, we turn an apparent limitation into 9 an opportunity by using the diversity of responses to short-lived climate forcers across the CMIP6 multi-model en10 semble to ensure our D&A is insensitive to structural uncertainty. Our framework is developed using a Pearl-causal 11 perspective, but forthcoming research now underway will apply the framework to in situ measurements using a 12 Granger-causal perspective. While the hypothesis-based framework and accompanying generalized D&A formula 13 we develop should be widely applicable, we include a strong caution that the hypothesis-guided simplification of 14 the formula for the historical climatic record of CONUS as described in this paper will likely fail to hold in other 15 geographic regions and under future warming. 16


Introduction
(DJF, MAM, JJA, and SON). Suppressed in the notation throughout is geographic location; this is (for now) generic 122 but will generally refer to geospatial locations (e.g., model grid cells or weather station locations). The formula is as including a sum of forced responses, anthropogenic or otherwise, and a set of cross-correlation terms between pairs 128 of forcings. (Note: the cross-correlation term is another way of writing a statistical interaction.) The forced term 129 P F (·) could also involve higher-order nonlinear terms (e.g., quadratic) or interactions (e.g., terms involving three 130 or more forcings). The terms F i (t, τ i ) are defined by Eq. A.4 (Appendix A) and represent forcings modified by the 131 lagged response of the climate system on characteristic timescales τ i determined by the thermal inertia of the ocean.

132
The internal variability component P I (·) is further decomposed as 133 P I (t) = P D (t) Low-freq. Drivers where P D (·) is a "driven" term that describes year-to-year variability due to known modes of large-scale ocean 134 (e.g., the El Niño Southern Oscillation or ENSO) or atmospheric (e.g., the Pacific North American pattern) drivers, 135 and P W (·) is a term associated with short-term weather variability. We assume that the driven component can be General Formula:  Finally, the internal variability component includes an error term P W (·) that represents residual variability in 144 the time series due to weather. We suppose that this term is, on average, zero, but has a variance described by 145 Var P W (t) = σ 2 (t) ≈ V 0 × exp{V 1 t}, i.e., the magnitude of the weather variability depends on the time t and can be approximately modeled as the product 146 of a "baseline" or pre-industrial variance V 0 and a time-varying quantity exp{V 1 t} where 1/|V 1 | is a characteristic t refer to annual quantities, we assume that all the temporal autocorrelation in the time series is fully captured 152 by P F (·) and P D (·), i.e., the residual weather variability terms for each year are temporally independent. For 153 the seasonal mean analysis, we assume that the weather variability follows a normal distribution (with mean zero 154 and variance σ 2 (t)), while for the seasonal Rx1Day analysis, we assume that the weather variability follows a 155 Generalized Extreme Value distribution (centered on zero such that the variability is described by σ 2 (t)). Note that 156 these assumptions on the distribution of weather variability refer to seasonal aggregates (mean and maximum) of States. Two notes: first, we set out to investigate a variety of complicated factors using the analyses described in 162 Section 4 and found that we could simplify Eq. 1 (the flow chart diagram in Figure 1 summarizes the various analyses 163 considered); the hypotheses are organized in terms of these discoveries. Second, it is important to emphasize that 164 our analyses focus on the contiguous United States over the historical record (1900 to present), hence all conclusions 165 regarding simplifications to the general formula apply only to a specific geographic region and time period. While 166 we do not claim that any of these hypotheses are satisfied in other parts of the globe or under future warming to 167 the global climate system, we assert that our hypothesis-driven framework and the accompanying general D&A 168 formula are broadly applicable to other regions and timeframes. To summarize the wide range of science questions 169 considered for each component of the general formula in Eq. 1, we have developed the flow chart shown in Figure 1. 170 Our first test addresses whether we can correctly identify the influence of WMGHG forcing on precipitation 171 using a linear regression framework: 172 Hypothesis 1 (H1). Globally, the coefficient for WMGHG forcing (β WMGHG ) in Eq. 2 is consistent with 2%/K for mean precipitation and 6%/K for extreme precipitation. 173 While the focus of this paper is CONUS, we do not have any a priori expectations about the magnitude of the 174 WMGHG signature regionally (e.g., over CONUS) and seasonally. However, we do know what to expect in terms 175 of global-average annual temperature changes, and hence our first test is conducted at the global scale. Next, again 176 looking at precipitation globally, we test whether there are any meaningful cross-correlation terms between external 177 forcings: 178 Hypothesis 2 (H2). Globally, over the historical period, all cross-correlation terms in the forced component in Eq. 2 are negligible, i.e., As mentioned previously, there could also be nonlinear components in the forced term given in Eq. 2; however, 180 for the current analysis we assume that these terms are negligible. The climate state over the historical period has 181 not changed enough to warrant strongly nonlinear behavior (e.g., step changes due to the loss of sea ice), and 182 numerous climate model studies have demonstrated the additivity of responses to forcing agents (e.g, Kirkevag et al.

183
[2008], Shiogama et al. [2013], Marvel et al. [2015]). Shifting our focus to CONUS, we next test for the presence of 184 meaningful trends in seasonal precipitation due to individual forcing agents: 185 Hypothesis 3 (H3). For CONUS and over the historical record, the effects of individual forcings are non-negligible for secular trends in precipitation. 186 In other words, this hypothesis seeks to determine the set of forcings that need to be included in the forced component 187 P F (·) of seasonal precipitation over CONUS and over the historical period. Given that one of these non-negligible 188 forcings is tropospheric aerosols (see Section 4.2), we next explore a set of tests to assess structural uncertainties in 189 our representation of anthropogenic aerosols: 190 Hypothesis 4a (H4a). For CONUS and over the historical record, the effect of SO 2 on seasonal precipitation dominates the influence of other anthropogenic aerosol species, namely black and organic carbonaceous aerosols and ammonia.
Hypothesis 4b (H4b). For CONUS and over the historical record, SO 2 emissions correlate with changes in precipitation as well as deposition of SO 2 by rainfall and column integrated SO 2 mass.
Hypothesis 4c (H4c). Asian and European aerosols do not affect US rainfall, and slow responses to aerosols (teleconnected from sea-surface temperatures (SSTs)) also do not affect precipitation.
Hypothesis 4d (H4d). For CONUS and over the historical record, a single CONUS-wide time series of SO 2 emissions yields the same result as using local estimates of SO 2 emissions.
We next test whether individual forcing agents influence the relationships between the modes of natural variability 192 and seasonal precipitation: 193 Hypothesis 5 (H5). For CONUS and over the historical record, all higher-order and interaction terms between the climate drivers and external forcing agents in Eq. 4 are negligible, i.e., ∂F j F j (t, τ j ) + · · · ≈ 0.

194
Because externally forced trends in seasonal precipitation over the historical record are small, we next test whether 195 we can distinguish the influence of WMGHGs from the influence of anthropogenic aerosols: 196 Hypothesis 6 (H6). For CONUS and over the historical record, the compensating errors in the effects of WMGHG and anthropogenic aerosol forcing are negligible, i.e., the estimated effects of WMGHGs and aerosols are (1) unbiased and (2) do not involve any compensating trade-offs.

197
Finally, we conduct a test of the influence of warming on the signal-to-noise ratio in seasonal precipitation:

198
Hypothesis 7 (H7). For CONUS and over the historical record, the signal-to-noise ratio of seasonal precipitation is a constant over time, i.e., Each of these hypotheses will be evaluated using a large set of climate model output (most from the CMIP6 multimodel ensemble) in a perfect data sense to test and guide fits that will eventually be applied to observations and also ensure our D&A is insensitive to structural uncertainty. The result of our explorations (described in the remainder of this paper) is that we can simplify the general formula defined by Eqs. 2-5 as follows: The time lag τ WMGHG is estimated in Appendix A. For reasons discussed in section 4.3.3, we assume that the response 200 of precipitation to SO 2 is dominated by fast processes and therefore set τ SO 2 ≈ 0.  Table 3 of Hodnebrog et al. [2013]. To characterize the WMGHG forcing overall, we will simply compute a total 270 WMGHG forcing (CO 2 +CH 4 +N 2 O+CFC-11+CFC-12), all in one value as a function of time; see Figure 2(a).

271
Since the response of precipitation to WMGHG forcings is dominated by the lagged response of SSTs to these  Anthropogenic aerosols. At the CONUS scale, we require a time series of anthropogenic aerosols that can be applied 282 to the observational record. A significant amount of the analysis in this paper specifically addresses this question (see 283 Section 4.3); the punchline is that for the CMIP6-historical, DAMIP hist-aer, and LUMIP hist-noLu experiments we 284 find that a single CONUS-wide time series of SO 2 emissions derived from observations sufficiently characterizes 285 the influence of anthropogenic aerosols on seasonal precipitation. Specifically, we utilize the CMIP6 forcing datasets

312
We now describe a set of analyses that test each of the hypotheses outlined in Section 1. Figure 1   Our first investigation involves hypotheses H1 and H2, which ask whether we can correctly identify the magnitude 316 of the WMGHG effect on precipitation via a linear regression framework (H1), and if so, can we further isolate 317 this WMGHG dependence in a noisy climate system with all forcings (H2)? The hist-GHG ensembles provide a 318 clean test of hypothesis H1, since the only external influence on precipitation in these simulations is the WMGHG 319 forcing. For H2, we subsequently use the historical runs to assess whether the WMGHG dependence is the same in 320 both scenarios or the WMGHG dependence is meaningfully influenced by other external forcings. For both of these 321 tests, we focus on the global average of annual mean and extreme precipitation; here, the global average for extreme 322 precipitation involves the grid cell average of annual maximum daily precipitation (Rx1Day).

323
To test H1, we conduct three regression analyses: first, we regress where Y hist-GHG m,e  level.

332
The estimated hist-GHG percent changes in precipitation (mean and return values) and scaling rates as a function 333 of the warming are shown in Figure 3 for each model along with 95% confidence intervals along both axes. The 334 expected 2%/K and 6%/K for means and extremes, respectively, are shown with black lines. Interestingly, the mean 335 scaling rates seem to be systematically less than 2%/K (although in some cases the confidence interval covers 2%/K); 336 the extreme scaling rate also is generally less than 6%/K, except for the IPSL runs. The scaling factors do not appear     this hypothesis, to be as general as possible, we simply assess changes in seasonal precipitation (mean and extreme) series overlap zero. It is also generally the case that the hist-noLu trends are non-distinguishable from the historical 414 trends. As such, we also conclude that LULCC does not significantly influence seasonal precipitation over CONUS.

415
In summary, we have found that WMGHG forcing tends to increase seasonal precipitation, anthropogenic 416 aerosol forcing tends to decrease seasonal precipitation, and the influence of natural forcings, stratospheric ozone, 417 and LULCC is minimal and can be, for our analysis, ignored. These results support our confidence assessments in    We now justify each of these conclusions in the next subsections using climate model output from the AerChemMIP 443 experiments described in Section 3.1.

457
The results of this analysis are summarized in Table 1

539
[2021], we simultaneously assess the influence of WMGHGs, SO 2 , and the set of five climate drivers discussed in 540 Section 3.3 as follows: for the mean (extreme) precipitation analysis, we allow the mean (location parameter) to      The tallied results are shown in for the bias metric, while the relative standard deviation tends to show slightly less agreement with hist-nat. The 599 agreement for pattern correlation is somewhat reduced, relative to the bias and relative standard deviation, although 600 in almost all cases the agreement is better than 50% (and in many cases much higher). This could be due in part    Once again, the single-forcing DAMIP experiments (hist-GHG and hist-aer) together with the CMIP6 historical 623 runs provide a test bed for answering this question. Again using the gridcell-specific analyses for these three ex-624 periments outlined in Section 4.4, note that across experiments we are modeling both x(t) and y(t) as a constant 625 multiplied by a time series of forcing (in the terminology of Eq. 6, this is x(t) = β WMGHG F WMGHG (t, τ WMGHG,m ) and 626 y(t) = β SO 2 F SO 2 (t, τ SO 2 )). Since the forcing time series are the same for the single forcing and historical runs, then if unbiased. Figure 7 shows the CONUS-average bias for the WMGHG and SO 2 coefficients for the single-forcing ex-630 periments relative to the historical (i.e., single forcing minus historical) along with 95% basic bootstrap confidence 631 intervals. For SO 2 , the bias is indistinguishable from zero for almost every model-season; for WMGHGs, the bias 632 is somewhat more distinguishable from zero (in that two to five of the models have confidence intervals that do not 633 overlap zero), but in most cases the confidence intervals include zero. For WMGHGs, when the bias is nonzero, it is 634 generally positive, meaning that we (on average) underestimate the influence of WMGHGs in the historical record.

635
Nonetheless, given this general agreement between the single-forcing and historical experiments, we can indeed 636 distinguish the influence of WMGHGs from the influence of anthropogenic aerosols.

637
The same coefficient estimates can be further interrogated to assess whether the single-forcing vs. historical 638 biases for WMGHGs are correlated with the SO 2 biases. Anti-correlation in the grid-cell WMGHG and SO 2 bias 639 would imply that there is a compensating trade off between our estimates of the effect of these two forcings. We 640 seek to minimize equal and opposite errors ǫ(t) in x(t) = β WMGHG F WMGHG (t, τ WMGHG ) and y(t) = β SO 2 F SO 2 (t, τ SO 2 )) 641 such that the fit to the observed precipitation time series is unperturbed, yet yields biased estimates of β WMGHG and 642 β SO 2 , i.e., 643 To answer this question, we calculate the Pearson correlation of the grid-cell SO 2 and WMGHG biases as well as    The signal-to-noise ratio of seasonal precipitation is a constant with respect to warming.

Yes
Virtually certain † Note: "emiso2" refers to SO2 emissions; "wetso2" refers to deposition of SO2 by rainfall; "iso2" refers to column integrated SO2 mass. In this paper, we have described novel methodology to address the difficult problem of conducting regional detection 689 and attribution for mean and extreme precipitation. Our framework can simultaneously account for anthropogenic 690 forcing (the "signal") while explicitly modeling the natural variability (the "noise"). Furthermore, we use output from 691 global climate models in a perfect data sense to rigorously test the application of a general formula for modeling 692 a time series of precipitation. Modeling uncertainty is well-known to be one of the main obstacles to conclusive 693 attribution statements for regional precipitation, but here we turn this limitation into an opportunity by making sure 694 the D&A methodology developed in this paper is insensitive to model uncertainty. Our focus in this paper has 695 been on the CONUS due to the relatively low observational uncertainty, with a forthcoming paper applying our 696 framework to century-length records of in situ measurements of seasonal precipitation. However, in principle, the series of hypotheses enumerated in this paper could be similarly applied to other global land regions to develop an 698 appropriate formula for regional D&A of precipitation.

699
The analyses described in this paper are extensive and involve a large set of global climate model output. We 700 have already attempted to summarize the scientific process that led to the development and testing of each hypothesis 701 with the flow chart shown in Figure 1. To further summarize our conclusions, we developed Table 3, which restates 702 the hypotheses, the model data sets used for each hypothesis, our conclusions, and a measure of our confidence 703 in each conclusion. The confidence statement attempts to summarize the evidence for each hypothesis, albeit in a 704 highly aggregated manner; nonetheless, we intend that a high-level statement for each of the lengthy analyses in 705 the paper will provide a helpful and concise summary of our conclusions. Specific information on how these labels formula, let us hypothesize that the regional precipitation P (t) tracks changes in global mean SST δT (t) as a function of time t and that 1031 lags in the precipitation response are due to time lags in the response of the global temperature: R = Unitless ratio of regional-to-global precipitation γ = Global precipitation sensitivity: 2%/K for mean and 6%/K for extremes The time evolution of δT (t) is dictated by the first law of thermodynamics applied to the global ocean: Define a timescale τ = ρ cp h/α. Then the first law can be written as: The solution is: In the limit τ = 0, one can show that, as expected from Eq. A.2 when the system is in equilibrium, This follows from Eq.s A.3 and A.4 if we define an integration variable x = (t − t ′ )/τ so that Substituting Eq. A.3 into Eq. A.1 gives an expression for the β coefficients that appear, for example, in Eq. 2 for the forced component: where we have introduced the possibility that β might exhibit dependence on τ .   A.2 Decreasing lagged forcing with increasing lag timescales then term B is also negative indefinite. Therefore If P (t) is prescribed and the forcing obeys the assumptions used to obtain Eq. A.7, then Eqs. A.5 and A.7 imply that 1043 d dτ The positive trend in β(τ ) with respect to τ is observed the DAMIP WMGHG-only runs ( Figure A.1).
This is indeed the behavior that's observed in climate models for t 20 years in Gregory regression plots of N (t) versus δT (t) [Gregory characterizing the "fast" and "slow" evolution of N i (t) and the "break" between these two regimes as follows: where Θ is the Heaviside function and F 0 = 8.053 is determined from the formulae in Etminan et al. [2016]. 1057 We use a Bayesian least squares fit to estimate the timescales τ f,i , τ s,i , and τ b,i in Equation A.15 using the N i (t) time series 1058 from each GCM with abrupt4xCO2 ensembles in Table E  year in both the piControl and 1pctCO2 runs. Furthermore, we derive the time series of WMGHG forcing for the 1pctCO2 runs using the 1072 specified boundary conditions for the 1pctCO2 experiment and various formulae described in Section 3.2. We then estimate the changes As described in Section 5, we provided a summary of the extensive analysis in this paper in Table 3, which includes a measure of 1119 our confidence in the conclusion of each hypothesis. Any attempt to distill the results of each analysis might feel oversimplified since 1120 the various hypothesis analyze numerous models, four seasons, and two types of precipitation (mean and extreme  implies "about as likely as not," ≥ 0.1 and < 0.33 implies "unlikely," ≥ 0.01 and < 0.1 implies "very unlikely," and < 0.01 implies 1132 "exceptionally unlikely." We now describe the specific quantities used to determine the confidence label for each hypothesis.   Here, "does not differ" is determined when the 95% confidence interval includes 0. As shown in Table 1  Given that both quantities result in a label of "likely," this is the category used in Table 3 for H6.