Archive for the ‘ecosystem model’ Category

Trees, volcanoes and climate

..large eruptions in the tropics and high latitudes were primary drivers of interannual-to-decadal temperature variability in the Northern Hemisphere during the past 2,500 years. Overall, cooling was proportional to the magnitude of volcanic forcing and persisted for up to ten years after some of the largest eruptive episodes. Our revised timescale now firmly implicates volcanic eruptions as catalysts in the major century pandemics, famines, and socioeconomic disruptions in Eurasia and Mesoamerica..[1]


  • a large volcanic eruption injects tens of millions of tons of sulphur dioxide into the stratosphere, causing global or hemispheric dimming. The cooling effect may last for several years.
  • the box-and-whisker graph above shows distributions of the annual growth of trees at two subarctic (taiga) forest sites in Quebec and  Russia and at a subalpine forest in France, from 1300 to the present day. These sites are close to the tree line where temperature is likely to have been the limiting factor for growth. The graph suggests that volcanic cooling had a significant impact on forest growth.
  • tree ring chronologies for each site were created by detrending raw data from individual trees and averaging. While this procedure loses information about climatic shifts lasting longer than a few decades, shorter term variability (for example, from volcanoes) is retained.
  • sulfur deposition in Antarctic ice cores was taken as the measure of volcanic aerosols. A simple annual average from 21 carefully dated Antarctic ices cores was taken[1]. To allow for the fact that aerosol effects can last for more than one year, an exponentially weighted moving average with decay constant of 3 years was computed. Years in which this average exceeded a threshold (40ng/g) were deemed “volcanic”. There were 30 such years since 1300, an average of 4 per century. Interestingly there have been no such years since 1821.
  • the list of “volcanic” years according to the above criterion since 1300 are: 1347 1459 1460 1461 1462 1463 1464 1465 1466 1600 1601 1602 1603 1604 1643 1644 1696 1697 1698 1810 1811 1812 1813 1815 1816 1817 1818 1819 1820 1821. Most of these correspond to years following well-known historical eruptions such as Kuwae 1453, Huaynaputina 1600, Tambora 1815 and the mystery 1809 volcano.
  • strangely, the impact of volcanoes appears weaker if Greenland ice core data is used instead of Antarctic data. Furthermore, the relation between tree-ring and ice core data appears much weaker in the pre-1300 data. For example, there is little evidence of the very large 1257 Samalas eruption in the growth indices at the sites selected (below).





Calculation Details

Tree ring indices were calculated from raw research data (.rwl files) archived by NOAA paleoclimate. For example, the Quebec l1 data is available here:

Growth index chronologies were calculated in R using the dendrochronology dplR package. Simple negative exponential detrending was used, which attempts to capture mean biological growth rates over the life of a tree (typically conifer ~ 100 years). For example, the raw Quebec l1 site index chronology (Gennaretti et al ) was calculated as follows:

Antarctic sulfate ice core data described in Sigl et al are available in this spreadsheet.

[1] Timing and climate forcing of volcanic eruptions for the past 2,500 years. Sigl, M., Winstrup, M., McConnell, J. R., Welten, K. C., Plunkett, G., Ludlow, F., … Woodruff, T. E. (2015).  Nature, 523(7562), 543–549. DOI: 10.1038/nature14565

Large-scale vegetation-rainfall correlations in Africa

Weather and climatic variability affect vegetation on continental scales. Intuitively, healthy vegetation is greener. A quantitative index of “greenness” called Normalised Difference Vegetation Index (NDVI) was developed by NASA.  The index takes values between 0 and 1. NDVI is related to the fraction of photosynthetically active radiation absorbed by the vegetation canopy. It is a powerful tool in vegetation monitoring.

Quantifying the relationship between weather and vegetation on a regional scale is a statistical problem which requires long time-series of NDVI. GIMMS is a 25 year (1981-2006) NDVI dataset based on inputs from a sequence of NOAA weather satellites carrying AVHRR radiometers. The GIMMS team removed satellite drift, calibration effects etc.[1] Of course NDVI has a complex spatial variation which is related to variations in vegetation cover type, not to weather variations. To isolate the impact of weather variations, it makes sense to “standardize” NDVI. Vegetation Condition Index (VCI) is:

latex VCI=\frac{NDVI \; \:-\; \:NDVI_{min}}{NDVI_{max}-NDVI_{min}}

For example, NDVImin might refer to the minimum local May NDVI value observed during 1981-2006. A tropical rainforest and a semi-arid region have very different NDVI, but both may have VCI close to 1 during a favourable period.

The maps below show the correlation between VCI and 3-month Standardized Precipitation Index (SPI3) computed for Africa by calendar month. GIMMS 8km NDVI was upscaled to 200km to match the scale of gridded monthly climate datasets for precipitation (GPCP).[2] The red areas correspond to areas of high correlation between VCI and SPI3. Vegetation in these areas responsive to rainfall variations on a 3-month timescale. A striking feature  is the line high correlation in the Sahel region. Vegetation in the Sahel is highly sensitive to drought. Southern Africa is also highly rainfall sensitive during Dec-May. It is impressive that the combination of datasets of very different origins (GPCP and GIMMS) produces spatially and temporally consistent information.

Sensitivity of Vegetation to Precipitation 1981-2006


R code snippets

Most of the effort is in pre-processing the climate and NDVI data so that they share the same projection, grid, time coordinate etc. GDAL and the interp() function from the akima package were used for this purpose. Eventually we end up with SPI3 and VCI data represented by N × T matrices where N = Nx ×Ny is the dimension of the spatial maps and T is the number of observations. From there,  SPI3 and VCI maps were organized by calendar month,

vci.m <- lapply(seq(1:12),function(m) vci[,seq(from=m, to=T, by=12)])
spi3.m <- lapply(seq(1:12), function(m) spi3[,seq(from=m, to=T, by=12)])

Correlation maps between VCI and SPI3  for each calendar month are obtained from,

corr.maps <- lapply(seq(1:12),function(m) sapply(seq(1:N), function(i) cor( vci.m[[m]][i,],spi3.m[[m]][i,]) ))

Finally the above map table was produced using spplot() from the sp package with overlays from the maps package.

[1]  Tucker, C.J., J.E. Pinzon, and M.E. Brown (2004), Global Inventory Modeling and Mapping Studies, NA94apr15b.n11-VIg, 2.0, Global Land Cover Facility, University of Maryland, College Park, Maryland, 04/15/1994.

[2] SG: Adler, R.F., G.J. Huffman, A. Chang, R. Ferraro, P. Xie, J. Janowiak, B. Rudolf, U. Schneider, S. Curtis, D. Bolvin, A. Gruber, J. Susskind, P. Arkin, 2003: The Version 2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979-Present). J. Hydrometeor., 4,1147-1167.