3.2.2. Model 1: Annual Central Valley Runoff Reconstruction from Tree-Ring Data

Separate models were developed to reconstruct 8RI annual flows over the periods 903– 2008 (long record) and 1640–2001 (short record). The long record prioritizes reconstruction *length* by making use of a small set of long tree-ring chronologies. The short record prioritizes reconstruction *accuracy* by taking advantage of a larger number of chronologies that, while not of great age, yield improved reconstruction accuracy and several centuries extension of flow beyond the gage record. The procedure described below was repeated for each of the two model periods.

A two-stage reconstruction method, introduced for reconstruction of river basin precipitation [94], and later extended for reconstruction of streamflow [53,95] was modified for this study. The first stage is conversion, by regression, of each of the available *N* chronologies into a separate single-site reconstruction (SSR) of *yt*, the square-root-transformed annual 8RI flows. A square root transform was found adequate to correct problems with violation of assumptions about the regression residuals that occur when using the untransformed flows as the regression predictand in the reconstruction models. The second stage, called multi-site reconstruction (MSR) combines the signals from the individual SSRs to get the final single time series of reconstructed flows. The two stages of reconstruction are outlined below and are described in more detail in Supplementary Material B.

In the first stage, *yt* is regressed stepwise on a pool of potential predictors that includes a site chronology, *xt*, its square, *x*<sup>2</sup> *<sup>t</sup>* , and lags *t*−2 to *t* + 2 on those two variables [96]. Preliminary stepwise modeling using cross-validation [97] is first applied to identify the step, *m*, beyond which addition of another predictor fails to increase the validation skill as measured by the reduction of error statistic (RE) [98]. The final SSR regression model has only those predictors entered in the first *m* steps, and substitution of the full available length of tree-ring predictors into the fitted equation gives the SSR reconstruction. Calibration accuracy of the SSR model is summarized by regression *R*2. Significance of the regression model is assessed by *p*F, the *p*-value of the overall-F of regression, and calibration uncertainty is summarized by the standard error of the estimate, or root-mean-square error (RMSEc) of calibration [96]. Validation accuracy is summarized by the root-mean-square error of cross-validation (RMSEv) which is computed from the cross-validation errors (see Supplementary Material B).

The product of the first stage of reconstruction in this study is a separate SSR, *y*ˆ*t*,*i*, of transformed flow for each of the *i* = 1, ... , *N* tree-ring chronologies (*N* = 69) mapped in Figure 1. Those whose SSR calibration signal is not significant (*p*<sup>F</sup> > 0.05) or whose SSR has no skill of validation (RE ≤ 0) were eliminated from the study. Depending on the time coverage (varies over SSRs), the remaining SSRs could contribute in the second stage of reconstruction.

The second stage is a re-calibration of the arithmetic mean of a subset of the *N* individual SSRs with acceptably strong signal and common time coverage into a final reconstruction—long or short, depending on the particular subset. The regression model for the MSR is

$$y = a + b\mathbf{x}\_{\mathbf{t}} + \mathbf{e}\_{\mathbf{t}\_{\mathbf{t}}} \tag{1}$$

where *xt* is the arithmetic average of the SSRs covering either the long or short records, *yt* is the square-root-transformed 8RI WY flow, *et* is the error term, and {*a*, *b*} are regression coefficients. The arithmetic average of the SSRs was preferred as a predictor with the assumption that the flow signal is best represented by the common variation in the SSRs. Since, by definition, the SSRs have variance proportional to the strength of their flow signals, a simple arithmetic average emphasizes chronologies with a strong flow signal. Equation (1) applies to both the long and short reconstructions, but for each *xt* is an average over a different set of SSRs. The calibration period is defined by the overlap of 8RI flows (and *yt*) with the particular SSR subset: WYs 1906–2008 for the long record and WYs 1906–2001 for the short record. Calibration and validation accuracy for the MSR models were measured by the same statistics already described for the SSR models. Substitution of the full-length SSRs into the fitted equations gave long and short reconstructions covering 903–2008 and 1640–2001, respectively. The RMSEv of the model was used, along with the assumption that the reconstruction residuals are normally distributed, to place a 50% confidence interval on the annual reconstructed transformed flows, *y*ˆ*t*. As a final step, the MSR reconstructions were back-transformed to original flow units (BCM) before interpretation.

As additional validation, time series plots and the Spearman correlation coefficient [99] were used to check agreement of the long and short reconstructions of 8RI flows spanning WYs 1872–1900 and reported in Bulletin 5 [65]. This step serves as a completely independent verification, as the Bulletin 5 data precede the start of the period used for screening treering data and calibrating and cross-validating the reconstruction models. Supplementary Material C provides the statistics of the SSR models.

Low-frequency features (decadal and longer) are of interest in understanding longterm hydrologic patterns. Severity of droughts and wet periods is also summarized by simple moving averages of reconstructed flows. Consistency of with other work was checked by comparing reconstructions with the sum of separate reconstructions of annual Sacramento River and San Joaquin River flow [20].

3.2.3. Model 2: Pre-Development Outflow and Salinity Reconstruction

The annual 8RI time series tree-ring reconstructions from Model 1 were used as the basis for reconstructing annual Delta outflow volume and salinity in the estuary under pre-development conditions. The modeling logic employed for the long- and short-period reconstructions is presented as a flow chart in Figure 2.

**Figure 2.** Flow chart for pre-development (Model 2) and contemporary (Model 3) models. <sup>1</sup> Model 2 calibration based on pre-development daily outflow (DWR, 2016) and daily X2-outflow relationship; Model 3 calibration based on historical daily and contemporary daily X2-outflow relationship; <sup>2</sup> Model 2 calibration based on pre-development annual outflow and historical Eight River Index (WYs 1922–2014); Model 3 calibration based on historical annual outflow and historical Eight River Index (WYs 1912–2018); <sup>3</sup> Tree-Ring reconstruction from Model 1 for long record and short record through WY 1850.

Annual pre-development Delta outflow volume was estimated for both reconstruction periods assuming a power law relationship between Delta outflow volume and annual Central Valley runoff volume:

$$
\Delta \text{ellta Outflow} = \alpha\_1 \times 8 \text{RI}^{\alpha\_2} \tag{2}
$$

where outflow and runoff volumes are in units of BCM per year and α<sup>1</sup> and α<sup>2</sup> are fitting parameters determined through least squares analysis [25] simulated pre-development conditions in the Central Valley and Delta assuming a historical runoff pattern measured over the 93-year period spanning WYs 1922–2014. Annual Delta outflow volume from the CDWR simulation, along with annual Central Valley runoff (as measured by the 8RI over the same 93-year period), were used to calibrate Equation (2). Tree-ring reconstructions of Central Valley runoff (i.e., the 8RI) were used in conjunction with Equation (2) to estimate pre-development annual Delta outflow volumes through WY 1850.

Seasonal (February–June and July–October) average pre-development X2 positions were estimated for the reconstruction periods assuming power law relationships between X2 position and annual Delta outflow volume:

$$\lambda \mathbf{2} = \alpha\_3 \times \text{Delta} \, \text{Out} \, \text{Out} \, \text{flow}^{\infty 4} \tag{3}$$

where X2 position is in units of km from Golden Gate and α<sup>3</sup> and α<sup>4</sup> are fitting parameters determined through least squares analysis. The 93-year simulated pre-development Delta outflow time series, as described above, was utilized to calibrate Equation (3). Several modeling steps were followed to generate 93-year X2 calibration data sets. First, daily outflow from the aforementioned CDWR [25] simulation was transformed into a daily "antecedent outflow" time series to represent flow time-history in the estuary [14,84]. This transformed outflow time series was used to generate a daily X2 time series using a predevelopment flow–salinity relationship reported by [26]. Finally, this synthetic daily X2 time series was averaged to develop seasonal average X2 calibration time series. Time series of annual pre-development outflow volume, estimated from the tree-ring reconstructions (Model 1), were used in conjunction with Equation (3) to estimate seasonal average X2 positions through WY 1850.
