Global Temperature Response to Volcanic Activity

# Global Temperature Response to Volcanic Activity

---

## Motivation and background

$$
\Delta T(t)=\int\_{\infty}^{t}\phi(t-s)\[f\_K(s)\mathrm{d}s+\sigma \mathrm{d}B(s)\]
$$

_Linear and stationary response $\phi(t)$_

Note:

Let us start off with an equation describing temperature fluctuations in terms of an
integral over a response function G, multiplited with a known forcing and a noise term.
This assumes a linear and stationary response, and we want to test the hypothesis of a
linear response function G.

### Filtered Poisson Process

$$
\Delta T\_K(t)=\sum_{k=1}^K A\_k \phi\left( \frac{t-t\_k}{\tau\_\mathrm{d}} \right)
$$

$$
\downarrow
$$

$$
\Delta T\_K(t)=\[\phi\*f\_K\]\left(\frac{t}{\tau_\mathrm{d}}\right)
$$

Note:

If only volcanic forcing is considered, we may write the equation for temperature
fluctuations on a different form, as a sum over K temperature time series that are the
response to K volcanic events, each with its own amplitude and arrival time, and with a
common width.

This formalism allows us to write this up as a convolution between a general response
function and a forcing time series. The amplitudes and arrival times are now baked into
F, while the duration time is given by the shape of the response function.

Note:

21 events (K=21)

### Deconvolution algorithm

$$
\phi^{(n+1)}=\phi^{(n)}\frac{(T\_K-\langle T\_K\rangle)\*\hat{f}\_K+b}{\phi^{(n)}\*f\_K\*\hat{f}\_K+b}
$$

<a href="https://doi.org/" data-citation-key="@richardson1972">Richardson (1972)</a>,
<a href="https://doi.org/10.1086/111605" data-citation-key="@lucy1974">Lucy (1974)</a>,
<a href="https://doi.org/10.1088/0266-5611/26/2/025004" data-citation-key="@benvenuto2009">Benvenuto
et al. (2009)</a>

Notes:

An iterative algorithm where an initial guess of the response function is updated, with
the times series for temperature T and forcing F also on the right-hand-side. The term b
in both the numerator and the denominator ensures positive values.

The hat on the forcing indicate that is uses negative times, i.e., F_hat(t)=F(-t).

### Model setup

1. Synthetic eruptions are given to CESM2
2. Run using the BWma1850 compset

Notes:

What we want is to run CESM2 with custom synthetic volcanoes while the rest of the
climate system stay unforced by external factors.

Made a Python project to create the volcanic forcing files. Uses the
**createVolcEruptV3.ncl** script from
[svn.code.sf.net](http://svn.code.sf.net/p/codescripts/code/trunk/ncl/emission/) to
convert from the source files.

The middle atmosphere (MA) chemistry is capable of calculating the evolution of
stratospheric aerosol from SO2 emissions from large volcanic eruptions. WACCM with MA
can only run at nominal 2° resolution, while the more expensive TSMLT1 (troposphere,
stratosphere, mesosphere and lower thermosphere) can go down to nominal 1° resolution.

---

## Results

### Single-volcano events

| <div class="fragment grow" data-fragment-index="0"><div class="fragment shrink" data-fragment-index="1"><img src="https://raw.githubusercontent.com/engeir/hack-md-notes/4c76fa84d73699f3dd51cf9a8234d9142e54e9d1/assets/pic/volcano-ensemble-waveforms/medium-waveform.png" ></div></div> | <div class="fragment grow" data-fragment-index="1"><div class="fragment shrink" data-fragment-index="2"><img src="https://raw.githubusercontent.com/engeir/hack-md-notes/4c76fa84d73699f3dd51cf9a8234d9142e54e9d1/assets/pic/volcano-ensemble-waveforms/medium-plus-waveform.png" ></div></div> | <div class="fragment grow" data-fragment-index="2"><img src="https://raw.githubusercontent.com/engeir/hack-md-notes/4c76fa84d73699f3dd51cf9a8234d9142e54e9d1/assets/pic/volcano-ensemble-waveforms/strong-waveform.png" ></div> |
| :----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------: | :---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------: | :-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------: |
| $26\\,\mathrm{Tg}$  | $400\\,\mathrm{Tg}$  | $1629\\,\mathrm{Tg}$  |

Notes:

The smallest eruption is comparable in size to the Mt. Pinatubo eruption. The Young(est)
Toba Tuff eruption (∼74000 years ago) has an estimated sulphur dioxide amount of 5400
Tg, more than the strongest simulated here.

Ensemble medians of four simulations with eruptions spaced with three months, one in
each season.

Notes:

We here look at the waveform of the temperature, to see if a clear difference between
the shapes, and in particular the tails, can be found. The two strongest eruptions
produce temperature responses that are similar in shape, while the weakest has a sharper
peak. The tail is more similar two the other two, but this could be an effect of the
normalization we are doing.

Note:

If we instead normalize by dividing through with the maximum value of the time series,
we find that the initial rise across all three is similar, while the weakest move faster
to equilibrium than the two stronger. One possibility is that the two stronger reach
past a threshold the weakest never get to.

### Double-volcano events

_Do non-linear effects become important as an eruption occur when the temperature is in
a perturbed state?_

Notes:

We see that when we superpose the single-events using the control run mean temperature
as the baseline, we largely are able to explain the shape of the double-event. Could not
have expected any better match, but this is clearly not enough to conclude.

Works well as a test to check if it is worthwhile doing.

Problem: simulations are not long enough to have a periodic signal

### CESM LME

CESM1 Last Millennium Ensemble

<a href="https://doi.org/10.1175/BAMS-D-14-00233.1" data-citation-key="@ottobliesner2016">Otto-Bliesner
et al. (2016)</a>

Notes:

The volcanic forcing is here events with zero width, meaning an event consist of one
element, described fully by an arrival time and amplitude.

Notes:

The temperature is the output of simulations with the CESM1 climate model.

Notes:

The estimated response function has the expected shape with a short initial rise and a
longer decay lasting approximately a decade.

A good amount of noise is still present, and adding to this analysis is a more in depth
look at how the deconvolution algorithm handle this kind of noisy input data.

### Deconvolution Issues

Notes:

Two things to notice:

1. Different resolution, so we must either up sample or down sample.
2. The forcing has a decaying tail, the events are not just a delta-like pulse with an
   amplitude and arrival time.

Note:

The corresponding response function has strong seasonal dependence, marked with filled
in circles, and otherwise falls down to zero very quickly.

Note:

This first example is the simplest case. Forcing with no width is convolved with a
response function (original response) which gives us a temperature time series.
Deconvolving with the forcing and temperature time series gives back an almost perfect
estimated response function.

But, what if we moved all forcing events one element forward before using it in the
deconvolution?

Forcing A

Note:

From this ...

(I.e., we use this to convolve with a response function and create the temperature.)

Forcing B

Note:

... to this.

(And this when we deconvolve.)

Note:

When the forcing events have no width, reconstructing the response function is fine. We
get some spurious spikes, but the general shape is correct

Note:

However, if the forcing events have a decaying tail, reconstructing the response fails.
It falls to zero immediately, except from a spike on the first index (time = 0).

---

## Going Forward

### Paper 1

Linear response to volcanic forcing in different climate states in CESM2

Note:

This relates to the first part of the results, with CESM2 simulations, with the
possibility of running the same experiments in a different climate state. (E.g., after a
CO2 doubling/quadrupling.)

### Paper 2

Pulse estimation from simulation data with a deconvolution method

Note:

This relates to the deconvolution of a long data set, and the challenges you meet when
using the deconvolution algorithm with climate data.

### Paper 3

Assessment of non-linear responses in the LongRunMIP ensemble

- LongRunMIP: ensemble of $1000$–$6000$ year long simulations
- Temperature anomaly as $$
  \Delta T(t)=\int\_{\infty}^{t}\phi(t-s)\[f\_K(s)\mathrm{d}s+\sigma \mathrm{d}B(s)\]
  $$
  
- Parameter estimation within a Bayesian framework, using INLA\*

\*integrated nested Laplace approximations, introduced for climate data in <a
href="https://doi.org/10.5194/esd-11-329-2020"
data-citation-key="@myrvollnilsen2020">Myrvoll-Nilsen et al. (2020)</a>

Note:

The LongRunMIP is an ensemble consisting of climate simulations between 1000 and 6000
years in length. Within this ensemble we want to take the model introduced in the very
beginning for temperature fluctuations to estimate the response from a general forcing.

To do this, parameter estimation is done within a Bayesian framework by making use of
the INLA methodology, which was presented by Myrvoll-Nilsen et al. (2020), who also
demonstrated its use on climate data.

---

## References

- 
<div class="csl-entry" id="ref-benvenuto2009" role="doc-biblioentry">
  Benvenuto, F, R Zanella, L Zanni, and M Bertero. 2009.
  <span
    >“Nonnegative Least-Squares Image Deblurring: Improved Gradient Projection
    Approaches.”</span
  >
  <em>Inverse Problems</em> 26 (2): 025004.
  <a href="https://doi.org/10.1088/0266-5611/26/2/025004"
    >https://doi.org/10.1088/0266-5611/26/2/025004</a
  >.
</div>
-

<div class="csl-entry" id="ref-lucy1974" role="doc-biblioentry">
  Lucy, L. B. 1974.
  <span
    >“<span class="nocase"
      >An iterative technique for the rectification of observed
      distributions</span
    >.”</span
  >
  <em>The Astronomical Journal</em> 79 (6): 745.
  <a href="https://doi.org/10.1086/111605">https://doi.org/10.1086/111605</a>.
</div>
-

<div class="csl-entry" id="ref-myrvollnilsen2020" role="doc-biblioentry">
  Myrvoll-Nilsen, E., S. H. Sørbye, H.-B. Fredriksen, H. Rue, and M. Rypdal.
  2020.
  <span
    >“Statistical Estimation of Global Surface Temperature Response to Forcing
    Under the Assumption of Temporal Scaling.”</span
  >
  <em>Earth System Dynamics</em> 11 (2): 329–45.
  <a href="https://doi.org/10.5194/esd-11-329-2020"
    >https://doi.org/10.5194/esd-11-329-2020</a
  >.
</div>
--
-

  <div class="csl-entry" id="ref-ottobliesner2016" role="doc-biblioentry">
  Otto-Bliesner, Bette L., Esther C. Brady, John Fasullo, Alexandra Jahn, Laura Landrum, Samantha Stevenson, Nan Rosenbloom, Andrew Mai, and Gary Strand. 2016. <span>“Climate Variability and Change Since 850 <span>CE</span>: An Ensemble Approach with the Community Earth System Model.”</span> <em>Bulletin of the American Meteorological Society</em> 97 (5): 735–54. <a href="https://doi.org/10.1175/BAMS-D-14-00233.1">https://doi.org/10.1175/BAMS-D-14-00233.1</a>.
  </div>
-

  <div class="csl-entry" id="ref-richardson1972" role="doc-biblioentry">
  Richardson, W. H. 1972. <span>“<span class="nocase">Bayesian-Based Iterative Method of Image Restoration*</span>.”</span> <em>J. Opt. Soc. Am.</em> 62 (1): 55–59. <a href="http://www.osapublishing.org/abstract.cfm?URI=josa-62-1-55">http://www.osapublishing.org/abstract.cfm?URI=josa-62-1-55</a>.
  </div>