Computing and visualizing transient global warming levels

# Atlas notebooks

***

> This notebook reproduces and extends parts of the figures and products of the AR6-WGI Atlas. It is part of a notebook collection available at https://github.com/IPCC-WG1/Atlas for reproducibility and reusability purposes. This work is licensed under a [Creative Commons Attribution 4.0 International License](http://creativecommons.org/licenses/by/4.0).
>
> ![Creative Commons License >](https://i.creativecommons.org/l/by/4.0/88x31.png)

Computing and visualizing transient global warming levels#

30/6/2021

J. Bedia - Santander Meteorology Group. Dept. of Applied Mathematics and Computing Science. University of Cantabria, Santander, Spain.

In this notebook we reproduce some of the key steps for the determination of transient Global Warming Levels (GWLs). The computation of the period corresponding to a given GWL is carried out by means of a external function available in the Atlas repository, which is applied in this notebook to a single model. The example illustrates the whole process of data access, data selection and preparation, computation of GWL periods and visualization.

Package dependencies#

This notebook is based on R programming language and requires packages:

magrittr to pipe (%>%) sequences of data operations improving readability

library(magrittr)

GWL calculation#

The function (warming-levels/scripts/getGWL.R) performs the key task of GWL determination. It is an atomic function to compute the timing of a user-defined transient Global Warming Level, that uses as input a numeric vector of mean global annual temperature projections for a given GCM. This numeric vector must be named, using the corresponding years as names. Both pre-industrial baseline (default to 1850-1900) and transient period (usually up to year 2100) are concatenated in a single input vector.

The following example is intended to work with the global mean (land+sea) 2-m temperature data stored as .csv files in the datasets-aggregated-regionally directory.

The getGWL function is loaded by issuing:

source("../warming-levels/scripts/getGWL.R")

The script is available for download:

getGWL.R

Example for one single CMIP5 GCM#

Reading input data#

The transient GWLs of +1.5, +2, +3 and +4\(^\circ\)C are here calculated for one single CMIP5 GCM. The EC-EARTH r12i1p1 ensemble member simulations forced by the RCP 8.5 scenario are picked as an example. Their data are located by searching with an adequate pattern into the data directory:

target <- "EC-EARTH"
modelfiles <- list.files("../datasets-aggregated-regionally/data/CMIP5/CMIP5_tas_landsea", 
                         full.names = TRUE,
                         pattern = target) %>% print()

[1] "../datasets-aggregated-regionally/data/CMIP5/CMIP5_tas_landsea/CMIP5_EC-EARTH_r12i1p1_historical.csv"
[2] "../datasets-aggregated-regionally/data/CMIP5/CMIP5_tas_landsea/CMIP5_EC-EARTH_r12i1p1_rcp26.csv"     
[3] "../datasets-aggregated-regionally/data/CMIP5/CMIP5_tas_landsea/CMIP5_EC-EARTH_r12i1p1_rcp45.csv"     
[4] "../datasets-aggregated-regionally/data/CMIP5/CMIP5_tas_landsea/CMIP5_EC-EARTH_r12i1p1_rcp85.csv"     

The target GCM (and/or ensemble member, if it applies) can be directly chosen using the appropriate target pattern string, as in the example code above. Note that CMIP5 and CMIP6 datasets live in different directories, following the same reference syntax.

First, the historical experiment is read. It will provide the baseline for the GWL calculation.

hist.regions <- grep("historical", modelfiles, value = TRUE) %>%
                read.table(header = TRUE, sep = ",", dec = ".", comment.char = "#") 
str(hist.regions)

Show code cell output Hide code cell output

'data.frame':	1872 obs. of  60 variables:
 $ date : Factor w/ 1872 levels "1850-01","1850-02",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ GIC  : num  -23.4 -23.8 -23.1 -17.5 -10.8 ...
 $ NWN  : num  -19.77 -19.2 -20.02 -9.23 -2.17 ...
 $ NEN  : num  -27.88 -27.58 -21.15 -12.95 -4.83 ...
 $ WNA  : num  -3.5 -1.19 -1.94 2.88 9.27 ...
 $ CNA  : num  -3.256 -0.767 1.275 8.376 14.322 ...
 $ ENA  : num  2.38 1.86 5.24 8.99 13.58 ...
 $ NCA  : num  12.8 13.6 14.5 18.8 20.3 ...
 $ SCA  : num  22.6 22.6 23.8 24.1 23.4 ...
 $ CAR  : num  22.4 22.1 22.7 22.8 23.3 ...
 $ NWS  : num  19.5 20.1 19.8 19.6 18.7 ...
 $ NSA  : num  21.7 22.4 22.2 22.1 21.2 ...
 $ NES  : num  21.9 21.8 21.7 21.6 21.2 ...
 $ SAM  : num  20 20.1 19.4 19.3 18.7 ...
 $ SWS  : num  15 15.7 14.8 14.2 12 ...
 $ SES  : num  22.1 22 20.6 19.4 15.7 ...
 $ SSA  : num  12.24 13.27 10.96 9.89 7.44 ...
 $ NEU  : num  -0.885 -1.989 -0.688 3.378 5.005 ...
 $ WCE  : num  -1.787 -0.202 3.175 6.775 11.896 ...
 $ EEU  : num  -13.17 -8.55 -2.09 1.8 12.37 ...
 $ MED  : num  7.32 8.66 10.59 11.55 16.02 ...
 $ SAH  : num  13.9 17.2 21 23 24.9 ...
 $ WAF  : num  22.3 23.9 24.7 23.6 23 ...
 $ CAF  : num  20.8 21.9 21.7 21.4 20.9 ...
 $ NEAF : num  21.3 22.8 22.2 23.9 23.6 ...
 $ SEAF : num  21.3 21 20.5 21 20.8 ...
 $ WSAF : num  19.6 19.8 19.5 18.9 16.5 ...
 $ ESAF : num  20.2 20.5 19.9 18.9 17.6 ...
 $ MDG  : num  23.1 23.4 23.1 22.8 22.3 ...
 $ RAR  : num  -30 -36.73 -24.41 -13.99 -3.24 ...
 $ WSB  : num  -22.6 -16.54 -9.09 -2.48 11.03 ...
 $ ESB  : num  -29.85 -24.53 -15.45 -5.43 5.19 ...
 $ RFE  : num  -25.21 -21.46 -12.54 -5.66 1.13 ...
 $ WCA  : num  -1.907 -0.103 4.884 9.228 17.403 ...
 $ ECA  : num  -14.9 -10.94 -5.11 2.3 8.83 ...
 $ TIB  : num  -15.24 -12.138 -10.491 -3.001 0.421 ...
 $ EAS  : num  1.47 4.01 7.14 13.37 16.18 ...
 $ ARP  : num  14.2 17 19.2 24.2 27.1 ...
 $ SAS  : num  15.2 17.8 20.3 25.3 25.7 ...
 $ SEA  : num  23.5 23.6 23.8 23.9 24 ...
 $ NAU  : num  25.6 25.2 24.7 24.7 23.4 ...
 $ CAU  : num  28 27.8 25.8 19.5 16.7 ...
 $ EAU  : num  21.8 22.6 21.1 18.3 15.9 ...
 $ SAU  : num  18.7 18.9 18.2 15 13.7 ...
 $ NZ   : num  16.3 17 15.8 14.2 13.1 ...
 $ EAN  : num  -16.3 -20 -27.3 -30.1 -34.1 ...
 $ WAN  : num  -6.67 -8.04 -13.27 -18.48 -18.4 ...
 $ ARO  : num  -24.9 -34.1 -29.7 -19.6 -10.2 ...
 $ NPO  : num  15.8 15.3 15.6 16.4 17.6 ...
 $ EPO  : num  23.8 24 24 24.1 24 ...
 $ SPO  : num  19.5 20 19.8 19 17.8 ...
 $ NAO  : num  14.6 14.3 14.3 14.9 15.9 ...
 $ EAO  : num  24 24.2 24.4 24.3 23.8 ...
 $ SAO  : num  19.4 20.3 20.4 19.4 17.5 ...
 $ ARS  : num  23.1 23.1 23.5 24.3 25.6 ...
 $ BOB  : num  23.8 24.1 24.7 25.5 25.3 ...
 $ EIO  : num  23.9 23.9 24.4 24.7 24.7 ...
 $ SIO  : num  21.8 22.3 22.3 21.7 20.6 ...
 $ SOO  : num  7.4 8.16 7.55 6.05 4.33 ...
 $ world: num  10.6 10.9 11.6 12.6 13.4 ...

This file contains regionally-aggregated (land+sea) monthly temperature data for all the IPCC AR6 reference regions. From this file, we just retain the first (date) and last (world) columns, average by years and prepare a single named vector with the years:

yrs <- hist.regions[["date"]] %>% gsub("-.*", "", .)
hist <- tapply(hist.regions[["world"]], INDEX = yrs, FUN = "mean", na.rm = TRUE)
names(hist) <- unique(yrs)

We extracted the year as a string from the date column and computed the yearly average of the world averaged data. The last command sets the years as names for the numeric vector. We can now have a first look at the historical experiment temperature for this model.

plot(x = as.integer(names(hist)),
     y = hist, 
     ylab = "global tas (degC)", xlab = "year",
     type = "l", col = "blue", las = 1)
grid()
title(paste(target, "historical tas"))

../_images/00f662bce7bfb8dfd7ea3b30d699fa9c1ac5822b04e1ad743120886348d4ea35.png

The same procedure is next carried out with the RCP 8.5 experiment data, adapting the pattern used in the regular expression to retrieve the corresponding data and repeating the rest of the steps:

rcp85.regions <- grep("rcp85", modelfiles, value = TRUE) %>% 
                 read.table(header = TRUE, sep = ",", dec = ".", comment.char = "#")
yrs <- rcp85.regions[["date"]] %>% gsub("-.*", "", .)
rcp85 <- tapply(rcp85.regions[["world"]], INDEX = yrs, FUN = "mean", na.rm = TRUE)
names(rcp85) <- unique(yrs)

plot(x = as.integer(names(rcp85)),
     y = rcp85, 
     ylab = "global tas (degC)", xlab = "year",
     type = "l", col = "blue", las = 1)
grid()
title(paste(target, "RCP8.5 tas"))

../_images/bba4a74dcd62a459adea2ae2965425a3acbe1e2684a4e4d5616cb9abc29d559a.png

We can now also plot the concatenated historical + RCP85 time series. A vertical dashed line is added to differentiate both experiments.

tas <- append(hist, rcp85)
plot(x = as.integer(names(tas)),
     y = tas, 
     ylab = "global tas (degC)", xlab = "year",
     type = "l", col = "blue", las = 1)
grid()
abline(v = 2006, col = "red", lty = 2)
title(paste(target, "hist+RCP 8.5 tas"))

../_images/28986fa47598734af3d8cb949402aa1a14eaf7289a5af90f1a243217c0e86e41.png

Transient GWL calculation#

For the IPCC WGI Interactive Atlas, the time periods for which +1.5, +2, +3 and +4 degree transient Global Warming Levels (GWLs) are reached (with respect to pre-industrial 1850-1900 mean value) are computed for CMIP5 and CMIP6 data using a 20-year moving window.

The default output of getGWL corresponds to the central year (n) of the 20-year period when the warming is first reached (the GWL period is thus [n-9, n+10]). This function will return ‘NA’ when the GWL is not reached before (the central year) 2100. Note that the 20-year window width can be adjusted by the user to other arbitrary lengths. In that case, the interval will adapt accordingly.

The tas time series depicted in the previous time-series plot will serve as input for the getGWL function. In the next code chunk, the +1.5 degC GWL is calculated, considering the default baseline, pre-industrial period (1850-1900) and moving window width (20 years) used in the IPCC WG1 Atlas:

 gwl_1.5 <- getGWL(data = tas, 
                   base.period = c(1850, 1900),
                   proj.period = c(1971, 2100),
                   window = 20,
                   GWL = 1.5) %>% print()

In this particular example, the function returns an integer value (2018), corresponding to the central year (n) of the 20-year window for the selected GWL period.

gwl_1.5

In addition, the period intervals (calculated in this case as [n-9, n+10]) are indicated by the interval attribute (2009-2028) and can be retrieved by:

attr(gwl_1.5, "interval")

Therefore, the calculation of several GWLs for this GCM/experiment is straightforward with just looping over the GWL parameter:

gwls <- c(1.5, 2, 3, 4)
central.years <- c()
for (i in 1:length(gwls)) {
    central.years[i] <- getGWL(data = tas, 
                               base.period = c(1850, 1900),
                               proj.period = c(1971, 2100),
                               window = 20,
                               GWL = gwls[i]) 
}

Finally, these warming levels are plotted on top of the GCM temperature series:

plot(x = as.integer(names(tas)),
     y = tas, 
     ylab = "global tas (degC)", xlab = "year",
     type = "l", col = "blue", las = 1)
grid()
title(paste(target, "hist+RCP 8.5 tas"))

## Indicate GWLs as vertical lines
colors <- 1:length(gwls)
line.type <- "dashed"
abline(v = central.years, col = colors, lty = line.type)
legend("topleft",title = "GWLs RCP8.5",
       paste0("+", gwls, " degC (", central.years, ")"),
       lty = line.type, col = colors)

../_images/d06b9abd2592bacafa60ec48c59733c4c1e3538448f92b44d9738040b5fad121.png

A more detailed script gwl_time_series_plots.R for representing the GWL spread of CMIP5 and CMIP6 multimodel ensembles with a flexible parameter setting is provided in the warming-levels/scripts directory.

Session Information#

sessionInfo()