On the Value of Paleo-Reconstructed Data: Understanding Long Term Hydroclimatic Variability and Its Implications to Water Resources Management

background-size: contain

---

# On the Value of Paleo-Reconstructed Data:<br>Understanding Long-Term Hydroclimatic Variability and<br>Its Implications to Water Resources Management

*Presented by:*

**Nguyen Tan Thai Hung**]
]

**Prof. Brendan Buckley**]]
]

<div id="box1" class="streamflow">
  <p id="text1">
    streamflow
  </p>
</div>

---

# Instrumental streamflow data are not long enough

---

# Short records led to overallocation of water rights

]

<br>

**The Colorado River Compact (1922)**

Measured discharge: 20.2 km<sup>3</sup>/year

Allocation &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;: 18.6 km<sup>3</sup>/year

Long term average &nbsp;: 16.3&ndash; 17.6 km<sup>3</sup>/year

.center.firebrick[The Colorado no longer reaches the Pacific Ocean.]

.center.firebrick[The Compact's measurement period (1905&ndash;1922) was the wettest period in four centuries.]

???

Measured: 16.4 MAF

Allocation: 7.5 MAF upper, 7.5 MAF lower

Long-term: 13.2, 13.5, 14.3

1 km3 = 0.81 MAF

1 MAF = 1.23 km3

---

# Tree rings are proxies of past climate

.center.small[*Glyptostrobus pensilis* (swamp cypress) from Laos, image courtesy of Prof. Gretchen Coffman]

---

# Dendrohydrology: We can build models to relate<br> tree rings to soil moisture and streamflow

.pull-left[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-3-1.svg" width="504" style="display: block; margin: auto;" />
]

---

# Dendrohydrology: We can build models to relate<br> tree rings to soil moisture and streamflow

.pull-left[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-4-1.svg" width="504" style="display: block; margin: auto;" />
]

.pull-right.center[
<img src="data:image/png;base64,#figures/c2-mada.png" width="100%"/>

The Monsoon Asia Drought Atlas (MADA)

Cook *et al.* (2010)
]

---

# Dendrohydrology has a long history<br>with water management

Hardman and Reil (1936)

]

Schulman (1945)

]

Stockton and Jacoby (1976)

] 
]

---

.left-footnote[Meko, D. M., & Woodhouse, C. A. (2011). Application of Streamflow Reconstruction to Water Resources Management. In M. K. Hughes, T. W. Swetnam, & H. F. Diaz (Eds.), Dendroclimatology: Progress and Prospects]

.large[<br><span style="color:firebrick">**Applications in water management are still limited in scope and effectiveness.**</span>]

---

.absolute-center[
<svg height="360" width="360">
  <circle cx="180" cy="180" r="140" stroke="darkorange" stroke-width="5" fill="none"/>
</svg> 
]

10.1002/2017WR022114

R package: **`ldsr`**
]]]

10.1002/essoar.10504791.1

R package: **`mbr`**
]]]

10.1029/2020WR027883

R package: **`pprR`** (in progress)
]]]

Publications planned

R package: **`shy`** (in progress)
]]]

---

.slant-left[
<svg height="360" width="360">
  <circle cx="180" cy="180" r="140" stroke="darkorange" stroke-width="5" fill="none"/>
</svg> 
]

10.1002/2017WR022114

R package: **`ldsr`**
]]]

10.1002/essoar.10504791.1

R package: **`mbr`**
]]]

10.1029/2020WR027883

R package: **`pprR`** (in progress)
]]]

Publications planned

R package: **`shy`** (in progress)
]]]

.slant-right.steelblue[

Galelli, Nguyen, Turner, and Buckley (submitted)

(Chapter 6)
]

---

<br>

# Chapter 2

# Linear Dynamical System Reconstruction

]

]

---

# The Chao Phraya River Basin is a third of Thailand's area

.pull-right.center[

.center.small[Station P.1, Nawarat Bridge, Chiang Mai. Photo: Rachel Koh]
]

---

# Streamflow correlates well with the MADA

.big-left[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-5-1.png" width="648" />
]

Streamflow reconstructions from drought atlases

.small.grey[

Graham, N. E., & Hughes, M. K. (2007). Reconstructing the Mediaeval low stands of Mono Lake, Sierra Nevada, California, USA. The Holocene, 17(8), 1197–1210. DOI: 10.1177/0959683607085126

Adams, K. D., Negrini, R. M., Cook, E. R., & Rajagopal, S. (2015). Annually resolved late Holocene paleohydrology of the southern Sierra Nevada and Tulare Lake, California. Water Resources Research, 51(12), 9708–9724. DOI: 10.1002/2015WR017850

Ho, M., Lall, U., & Cook, E. R. (2016). Can a paleodrought record be used to reconstruct streamflow?: A case study for the Missouri River Basin. Water Resources Research, 52(7), 5195–5212. DOI: 10.1002/2015WR018444

Ho, M., Lall, U., Sun, X., & Cook, E. R. (2017). Multiscale temporal variability and regional patterns in 555 years of conterminous U.S. streamflow. Water Resources Research, 53(4), 3047–3066. DOI: 10.1002/2016WR019632
]
]

---

# Flood generation mechanism:<br> excess rainfall on wet soil

Stein, L., Pianosi, F., & Woods, R. (2020). Event‐based classification for global study of river flood generating processes. Hydrological Processes, 34(7), 1514–1529. DOI: 10.1002/hyp.13678

Lim, H. S., & Boochabun, K. (2012). Flood generation during the SW monsoon season in northern Thailand. Geological Society, London, Special Publications, 361(1), 7–20. DOI: 10.1144/SP361.3

]

---

# Most reconstructions do not account for catchment dynamics

.firebrick.large[**Conventional methods rely on linear regression**]

$$
y_t = \alpha + \beta u_t + \varepsilon_t
$$
]

$$
`\begin{equation}
  y\\
  u\\
  \alpha, \beta\\
  \varepsilon\\
\end{equation}`
$$
]

.right90.narrow50[

streamflow

principal components of MADA grid points

regression coefficients

white noise
] 
]

.pull-right.large.steelblue[

<br>

* How do we model catchment dynamics?

* What insights can we gain with a dynamic model?
]

---

# Linear Dynamical System (LDS)

---

# Linear Dynamical System (LDS)

---

# Linear Dynamical System (LDS)

.math[
$$
`\begin{align}
  x_{t+1} &= Ax_t + Bu_t + w_t \\
  y_t     &= Cx_t + Du_t + v_t \\
  w_t & \sim \mathcal{N}(0, Q) \\
  v_t & \sim \mathcal{N}(0, R) \\
  x_1 & \sim \mathcal{N}(\mu_1, V_1)
\end{align}`
$$
]

Collectively,

`$$\theta = (A, B, C, D, Q, R, \mu_1, V_1)$$`

]

<br>

$$
`\begin{align}
  x &\in \mathbb{R}^{p} \\
  y &\in \mathbb{R}^{q} \\
  u &\in \mathbb{R}^{m} \\
  A &\in \mathbb{R}^{p \times p} \\
  B &\in \mathbb{R}^{p \times m} \\
  C &\in \mathbb{R}^{p \times p} \\
  D &\in \mathbb{R}^{p \times p} \\
  Q &\in \mathbb{R}^{p \times p} \\
  R &\in \mathbb{R}^{p \times q} \\
\end{align}`
$$

]

.big-right.narrow60[
system state

system output (streamflow)

system input (PCs of MADA grid)

state transition matrix

input-state matrix

observation matrix

input-observation matrix

state-noise covariance

observation noise covariance  
]
]

---

# Learned with Expectation-Maximization

Need to handle missing data &rarr; .steelblue[Modified M-step (Lemma 2, page 20)]

.steelblue[
R package **`ldsr`** (available on CRAN and at [github.com/ntthung/ldsr](https://github.com/ntthung/ldsr))
]

---

# Catchment state improves flow estimation

.big-left[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-7-1.png" width="648" style="display: block; margin: auto;" />
]

<br>

$$
`\begin{align}
x_{t+1} &= Ax_t + Bu_t + w_t \\
y_t &= Cx_t + Du_t + v_t\\
C &> 0
\end{align}`
$$
]

]

<div id="htmlwidget-a9780af8fd51ab7ee65b" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-a9780af8fd51ab7ee65b">{"x":{"filter":"none","data":[["Linear regression","LDS"],[0.45,0.77],[0.41,0.4],[0.3,0.31],[0.48,0.53]],"container":"<table class=\"display\">\n  <thead>\n    <tr>\n      <th>Model<\/th>\n      <th>R²<\/th>\n      <th>RE<\/th>\n      <th>CE<\/th>\n      <th>KGE<\/th>\n    <\/tr>\n  <\/thead>\n<\/table>","options":{"dom":"t","columnDefs":[{"targets":1,"render":"function(data, type, row, meta) {\n    return type !== 'display' ? data : DTWidget.formatRound(data, 2, 3, \",\", \".\");\n  }"},{"targets":2,"render":"function(data, type, row, meta) {\n    return type !== 'display' ? data : DTWidget.formatRound(data, 2, 3, \",\", \".\");\n  }"},{"targets":3,"render":"function(data, type, row, meta) {\n    return type !== 'display' ? data : DTWidget.formatRound(data, 2, 3, \",\", \".\");\n  }"},{"targets":4,"render":"function(data, type, row, meta) {\n    return type !== 'display' ? data : DTWidget.formatRound(data, 2, 3, \",\", \".\");\n  }"},{"className":"dt-right","targets":[1,2,3,4]}],"order":[],"autoWidth":false,"orderClasses":false}},"evals":["options.columnDefs.0.render","options.columnDefs.1.render","options.columnDefs.2.render","options.columnDefs.3.render"],"jsHooks":[]}</script>

$$
`\begin{equation}
RE = 1 - \frac{\displaystyle\sum_{t \in \mathcal{V}} (y_t - \hat{y}_t)^2}{\displaystyle\sum_{t \in \mathcal{V}} (y_t - \bar{y}_c)^2}
\end{equation}`
$$
]

.math[
$$
`\begin{equation}
CE = 1 - \frac{\displaystyle\sum_{t \in \mathcal{V}} (y_t - \hat{y}_t)^2}{\displaystyle\sum_{t \in \mathcal{V}} (y_t - \bar{y}_v)^2}
\end{equation}`
$$

]

<br>

.math[
$$
`\begin{equation}
    KGE = 1 - \sqrt{(\rho - 1)^2 + \left(\frac{\hat{\mu}}{\mu} - 1\right)^2 + \left(\frac{\hat{\sigma}}{\sigma} - 1\right)^2}
\end{equation}`
$$
]
]

] 
]

---

# LDS reveals a history of regime shifts in the Ping River Basin

---

# Conclusions

* Streamflow reconstruction based on linear dynamical system

* Learned with modified Expectation-Maximization algorithm

]

* LDS reveals a history of regime shifts in the Ping River

* It is important to model catchment dynamics
]
]

---

<br>

# Chapter 4

# Monsoon Asia:

# Eight Centuries of Streamflow History

]

---

# MADA v2 is a significant upgrade from v1

.pull-left[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-11-1.png" width="504" />
]

<br>

<div id="htmlwidget-dafc9d3ece35786c38ea" style="width:100%;height:auto;" class="datatables html-widget"></div>
<script type="application/json" data-for="htmlwidget-dafc9d3ece35786c38ea">{"x":{"filter":"none","data":[["Resolution","Grid points","Chronologies","Time span"],["2.5° x 2.5°","534","327","1300–2005"],["1.0° x 1.0°","2724","453","0–2012"]],"container":"<table class=\"display\">\n  <thead>\n    <tr>\n      <th> <\/th>\n      <th>MADA v1<\/th>\n      <th>MADA v2<\/th>\n    <\/tr>\n  <\/thead>\n<\/table>","options":{"columnDefs":[{"className":"dt-center","targets":[1,2]},{"orderable":false,"targets":0}],"dom":"t","order":[],"autoWidth":false,"orderClasses":false}},"evals":[],"jsHooks":[]}</script>
]

---

# Grid points are selected by hydroclimatic similarity

.left-code[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-13-1.png" width="360" />

.small.grey[Knoben, W. J. M., Woods, R. A., & Freer, J. E. (2018). A Quantitative Hydrological Climate Classification Evaluated with Independent Streamflow Data. *Water Resources Research*, 54(7), 5088–5109. DOI: 10.1029/2018WR022913.]
]

<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-15-1.png" width="540" style="display: block; margin: auto;" />
]

---

# Reconstructions are reliable

.panelset[
.panel[.panel-name[Good skills at most stations]
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-16-1.png" width="936" style="display: block; margin: auto;" />

]

]
]

---

# Reconstructions capture historical events<br>and reveals spatial coherence of streamflow

.mid-right2.narrow60[
(1) Samalas eruption

(2) Angkor Drought I

(3) Angkor Drought II

(4) Kuwae eruption

(5) Ming Dynasty Drought

(6) Strange Parallels Drought

(7) East India Drought

(8) Tambora eruption

(9) Victorian Great Drought
]

---

# Asian rivers share common oceanic teleconnections

.small.bottom-right[Data for 1855&ndash;2012 <br> <br>]

---

# Teleconnections changed through time

.tall-fig[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/cor3-1.png" width="828" style="display: block; margin: auto;" />
]

---

# Conclusions

* Climate-informed automated grid point selection

]

* Streamflow in Monsoon Asia is spatially coherent

- .darkorange[Inter-basin transfer of water-dependent resources may have inadvertent effects]

* Common oceanic teleconnections

* Teleconnection changed through time

- .darkorange[Opportunity and challenge for climate-informed operations and forecasts]

]

---

<br>

# Chapter 3

# Multi-Proxy, Multi-Season Streamflow Reconstruction

]

]

---

# There have been few attempts on sub-annual reconstructions

* Stochastically disaggregate from annual reconstructions (Prairie *et al.*, 2008; Sauchyn and Ilich, 2017)

* With regional proxies as additional predictors (Stagge *et al.*, 2018)

<span style="color:firebrick">**Limitations**</span>

* Assume constant relationship between subannual and annual time series

* Do not fully leverage multi-sensory signals of proxies

] 
]

.center.huge[
<br>
Leverage proxy diversity

*directly* reconstruct subannual flow
]

<br>

.center.large[
<span style="color:steelblue">e.g., dry season, wet season, and water year reconstructions.</span>

] 
] 
]

---

# Challenges

1. We need a rich proxy network

2. How do we combine proxies?

3. How do we account for mass balance?
]

&nbsp; &nbsp; &nbsp; &nbsp; Total seasonal flow matches the annual flow closely

---

# The Southeast Asian Dendrochronology Network

<br>

**Tree ring cellulose stable oxygen isotope ratio**

`$$r_{sample} = \frac{^{18}\text{O}}{^{16}\text{O}}$$`

`$$\delta^{18}\text{O} = \left(\frac{r_{sample}}{r_{standard}} - 1\right) \times 1000$$`

]

---

# Annual mass balance

---

# Mass-balance-adjusted regression

Say we have the following regression equation:

`$$\mathbf{Y} = \mathbf{X}\beta + \varepsilon$$`

In reconstructions, `$\mathbf{Y}$` is streamflow and `$\mathbf{X}$` is a matrix of proxies.

We can form the least squares objective function

`$$\min_{\beta} \quad (\mathbf{Y - X}\beta)'(\mathbf{Y - X}\beta)$$`

and solve for `$\beta$` analytically:

`$$\beta = (\mathbf{X'X})^{-1}\mathbf{X'Y}$$`

]

.math[
$$
`\begin{align}
  \mathbf{Y_D} &= \mathbf{X_D\beta_D} + \mathbf{\varepsilon_D} \\
  \mathbf{Y_W} &= \mathbf{X_W\beta_W} + \mathbf{\varepsilon_W} \\
  \mathbf{Y_Q} &= \mathbf{X_Q\beta_Q} + \mathbf{\varepsilon_Q} \\
\end{align}`
$$
]

Let

.math[
$$
`\begin{align}
  \mathbf{Y} = \left[\begin{matrix}
        \mathbf{Y_D}\\
        \mathbf{Y_W}\\
        \mathbf{Y_Q}
      \end{matrix}\right], \quad
  \mathbf{X} = \left[\begin{matrix}
        \mathbf{X_D} &              & \\
                     & \mathbf{X_W} & \\
                     &              & \mathbf{X_Q}
      \end{matrix}\right], \quad
  \beta = \left[\begin{matrix}
              \mathbf{\beta_D}\\
              \mathbf{\beta_W}\\
              \mathbf{\beta_Q}
            \end{matrix}\right], \quad
  \varepsilon = \left[\begin{matrix}
        \mathbf{\varepsilon_D}\\
        \mathbf{\varepsilon_W}\\
        \mathbf{\varepsilon_Q}
      \end{matrix}\right]
\end{align}`
$$
]

We get back the canonical form:

It is tempting to impose a mass balance constraint:

`$$\mathbf{X_D\beta_D} + \mathbf{X_W\beta_W} = \mathbf{X_Q\beta_Q}$$`

$$
`\begin{align}
\delta &= \mathbf{X_D\beta_D} + \mathbf{X_W\beta_W} - \mathbf{X_Q\beta_Q} \\
       &= \mathbf{A}\beta
\end{align}`
$$

where `$\mathbf{A} = \left[\begin{matrix} \mathbf{X_D} & \mathbf{X_W} & -\mathbf{X_Q}\end{matrix}\right]$`

]

We introduce a penalty term with a weight `$\lambda$`

.math[
$$
\min_{\beta} \quad J = (\mathbf{Y - X}\beta)'(\mathbf{Y - X}\beta) + \lambda (\mathbf{A}\beta)'\mathbf{A}\beta
$$
]

The analytical solution is

] 
]

---

# Proxy combination

---

# Optimal proxy combination

Represent proxy combination with a binary vector

.math[
$$
`\begin{align}
  \mathbf{p} &= [p_1, p_2, ..., p_K]'\\
              & \\
  p_i &= \begin{cases}
            1 \qquad \text{if the } i^{\text{th}}  \text{ proxy is selected} \\
            0 \qquad \text{otherwise}
          \end{cases}
\end{align}`
$$
]

<br>

Each `$\mathbf{p}$` has a penalized least squares value `$J$`. We find `$\mathbf{p}$` with the smallest `$J$`.

This can be solved with *Genetic Algorithm*.

---

# The two models produce similar reconstructions...

---

# ...but Model 1 preserves the annual mass balance better

???

For Model 0, 28% of the time, the mass difference is outside the +- 10% MAF range, while for Model 1, that number is only 13%.

For Model 0, the difference can be as large as 640 Mm3, which is 90% of the Irrigation Demand of the Ping River downstream of Bhumibol.

---

# Conclusions

* Subannual reconstructions

* Optimal proxy combinations

* Mass balance adjusted regression

]

* Other climate targets, e.g. rainfall

* Tributaries' flows add up to main stream's flow

* Probabilistic assessment of water resource systems

* Higher resolutions (e.g., monthly)

]]

---

<br>

# Chapter 5

# Probabilistic assessment of water system <br> with monthly reconstructions
]

<img src="data:image/png;base64,#figures/sirikit-dam-cropped.png" width="100%"/>
]

---

# Chao Phraya water system simulation

**Workflow**

1. Reconstruct monthly inflow at P.1 and N.1

1. Create bias-corrected reconstructions with quantile mapping

1. Generate three stochastic ensembles: instrumental, reconstruction, bias-corrected reconstruction

1. Convert to inflow

1. Simulate

1. Assess the outputs

]

---

# From upstream reconstruction to stochastic inflow

**Step 2:** Bias correction with quantile mapping

$$
U_{bc} = F_o^{-1}(F_r(U_r))
$$

**Step 3:** Generate stochastic ensembles of upstreamflow

Vector AutoRegressive Moving Average (VARMA)

$$
`\begin{equation}
 \mathbf{y}_t = \sum_{i=1}^p \phi_p \mathbf{y}_{t-i} + \sum_{j=1}^q \theta_j \varepsilon_{t-j} + \varepsilon_t
\end{equation}`
$$

Model order determined with Bayesian Information Criterion: `$p = 1, q = 1$`.

**Step 4:** transform to inflow

$$
`\begin{equation}
Q^k_{n, m} = a^k_m U^k_{n, m} + b^k_m
\end{equation}`
$$

---

# Assessment criteria (step 6)

`$$H^k = \frac{1}{N} \sum_{n=1}^N \left( \sum_{m=1}^{12} \eta \rho g h(S^k_{n,m}) R^k_{n,m}\right)$$`

$$
`\begin{equation}
  W = \frac{1}{N}\sum_{n=1}^N \left(
    \frac{1}{12}\sum_{m=1}^{12} \min \left\{ 
      \frac{r^{Bhumibol}_{n,m} + r^{Sirikit}_{n,m}}{D_m}, 1
    \right\}
  \right)
\end{equation}`
$$

$$
`\begin{equation}
    N_s = \sum_{n=1}^N \sum_{m=1}^{12} \mathbb{1}_{\{s_{n,m} > 0\}}
\end{equation}`
$$

---

# Monthly upstream flow reconstructions

]

<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-23-1.png" width="936" style="display: block; margin: auto;" />
]
]

---

# Reconstruction reveals rich "internal" variability

.mid-right.firebrick[To bias-correct or not, that is the question.]

---

# Three ensembles produce three distributions of<br>system performance

.tall-fig2[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-26-1.png" width="648" style="display: block; margin: auto;" />
]

---

# How do we reconcile the three ensembles?

- Non-linear methods, e.g., LDS

- Add more proxies (blue light reflectance, wood density, etc.)
]

.left-footnote[Patskoski, J., & Sankarasubramanian, A. (2015). Improved reservoir sizing utilizing observed and reconstructed streamflows within a Bayesian combination framework. Water Resources Research, 51(7), 5677–5697. DOI: 10.1002/2014WR016189
]

.large.steelblue[**It is important to incorporate paleo data into water system assessment.**]

---

# Synthesis

.large.steelblue[
**Some important gaps are addressed**

* Catchment dynamics

* Monthly reconstructions

* Proxy selection methods

* Large scale variability in Monsoon Asia
]
]

* Reconciling distributional differences

* .darkorange[Climate change and variability in the future]
]
]

.large.center[Is it time to apply streamflow reconstructions to water management?]

.large.center[Yes!]

.large.pull-down-center[Thank you!]

---

# Extra slides

---

# Simultaneous Learning-Reconstruction

.big-fig[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-27-1.png" width="1152" />
]

.pull-down.large[Need to handle missing data &rarr; .steelblue[Modified M-step (Lemma 2, page 20)]]

---

# LDS serves as a stochastic streamflow generator

---

# Ensemble point-by-point regression with cross-validation

Thirty models are considered:

* .steelblue[Five climate distances

$$
d = 0.10, 0.15, 0.20, 0.25, 0.30
$$

]

* Six PCA weights (Cook *et al.*, 2010)

$$
z_i = g_i \rho_i^p
$$

$$
p = 0, \frac{1}{2}, \frac{2}{3}, 1, \frac{3}{2}, 2
$$
]

* hold-out 25%

* contiguous chunks

* repeat 50 times

* (CE, KGE) nearest to (1, 1)
]
]

<br>

---

# Counter example

---

# MBR full reconstruction

---

# Reconstructions reveal intra-annual variability

---

# Storage follows range of variability

---

# Hypothetical reservoir operating rules (steps 5 and 6)

.pull-left[
<img src="data:image/png;base64,#thesis-presentation_files/figure-html/unnamed-chunk-33-1.png" width="576" />
]

]