In brasmus/preciparea: Precipitation area ananlysis.
knitr::opts_chunk$set(echo = TRUE)
Statistics on the precipitation area
The instantaneous area of precipitation is an important indicator describing the hydrological cycle and the state of earth's climate. This is a document of data analysis carried out to elucidate this aspect of the climate, based on model data and historical observations, and is intended as the supporting material for the paper with the title 'Was the 1998--2016 decrease in the tropical precipitation area real?'. It is meant to provide a record like a lab note-book, to enhance the possibility of replication and provide details about the analysis. It provides openness and transparency, and is not mean to be a 'well-structured' document.
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
Getting the data
The TRMM data was downloaded (February 2017) from NASA data portals https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/ using he instructions provided in https://disc.sci.gsfc.nasa.gov/recipes/?q=recipes/How-to-Download-Data-Files-from-HTTP-Service-with-wget and https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/doc/TRMM_Readme_v3.pdf. It is necessary to create a user account and set up user details/coockies as descibed in the intructions. The code presented here was written for Linux, and it is uncertain whether it works for Windows machines. It maybe will work on iOS. The data volumes are large and require some time and space for a complete download
Satellite data - the Tropical Rain Measurement Mission (TRMM)
Download the daily TRMM data (~16Gb and takes time) on a Linux platform (using a terminal and bash-scripts).
```{bash, eval=FALSE}
Retrieving the data:
https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/
See this web page for how to download the files.
https://disc.sci.gsfc.nasa.gov/recipes/?q=recipes/How-to-Download-Data-Files-from-HTTP-Service-with-wget
https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/doc/TRMM_Readme_v3.pdf
Script for downloading thr TRMM data
!/bin/bash
mkdir TRMM
cd TRMM
for year in {1998..2016}
do
for month in {1..12}
do
wget --load-cookies ~/.urs_cookies --save-cookies ~/.urs_cookies --keep-session-cookies
-r -c -nH -nd -np -A nc4
"https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/$year/$month"
done
done
#### Independent reanalysis - ERAINT
The ERAINT data were downloaded from <http://apps.ecmwf.int/datasets/data/interim-full-moda/levtype=sfc/> (Monthly means of Daily means and daily precipitation) as netCDF and renamed `ERAINT_t2m.nc` (2-meter temperature), `ERAINT_Qs.nc` (surface moisture flux), and `ERAINT_Q-columntotal.nc` (total column water). Te daily precipitation was downloaded in multiple files due to large data volumes: `ERAINT_24hr_precip*.nc`. The unit of the precipitation was 'm' (it is assumed it was m/day).
The following lines of code were the contents of the script `totprecip.jnl` implemnted by PMEL Ferret:
```{bash, eval=FALSE}
use $1
set reg/y=50S:50
CANCEL LIST/HEAD
let TotP = TP[x=@din,y=@din]
let pre001 IF TP GT 0.0001 then 1 else 0
let Ap001 = pre001[x=@din,y=@din]
let pre01 IF TP GT 0.001 then 1 else 0
let Ap01 = pre01[x=@din,y=@din]
let pre1 IF TP GT 0.01 then 1 else 0
let Ap1 = pre1[x=@din,y=@din]
let pre0 IF TP LT -9999 then 0 else 1
let TotA=pre0[x=@din,y=@din]
!save/append/clobber/file=totprecip.nc TotP,Ap01,Ap1,TotA
list/append/file=totprecip.txt TotP, Ap001,Ap01,Ap1,TotA
cancel data ($1)
cancel memory
! Earth's area = 5.111859e+14 m^2
The scripty was ran through a bash job:
```{bash, eval=FALSE}
!/bin/bash
for f in ERAINT_24hr-precip*.nc; do
echo $f
/home/rasmusb/ferret/bin/ferret -gif -script ~/Downloads/totprecip.jnl $f
done
The use of MERRA-2 was also attempted, but without success. The right data (daily precipitation) was hard to find, and a search on the web site gave numerous hits with all sorts of precipitation results but not for 24-hr precipitation on a daily basis through the reanalysis interval. A Google search did not give any useful tips. An email request was never answered.
### The R-code
To keep all R-scripts/functions, extracted data and documentation together, a dedicated R-package `preciparea` was created that must be installed and activated for replicating this work. It is available at <https://github.com/brasmus/preciparea>.The R-scripts are open-code in this package and is part of the supporting material (SM) and contains both R-scripts and extracted data for this analysis.
The R-package can be installed from github:
```r
install.esd <- ("esd" %in% rownames(installed.packages()) == FALSE)
install.preciparea <- ("preciparea" %in% rownames(installed.packages()) == FALSE)
install.devtools <- ("devtools" %in% rownames(installed.packages()) == FALSE) &
(install.esd | install.preciparea)
if (install.devtools) {
print('Need to install the devtools package')
## You need online access.
install.packages('devtools')
}
library(devtools)
if (install.esd) install_github('metno/esd')
if (install.preciparea) install_github('brasmus/preciparea')
The package must be activated in R:
library(preciparea)
library(esd)
sessionInfo()$otherPkgs
Preparing indices from the data
This section shows how indeces stored in the preciparea package were derived, such as data("Parea").
TRMM
The TRMM netCDF files were unfortunately not stored with CF convention and a time dimension, and it was unclear how common netCDf processing tools like ncra could be applied to estimate the monthly means. Instead, the monthly data were aggrigated from the daily ones in the R-environment and the R-package preciparea (which was time-consuming).
The following command line were used to prepare the daily indicex describing the precipitation area $A_P = \int_A H(x-x_0)da$, where $H$ is the Heaviside function, $x$ is the precipitation, and $x_0$ is a precipitation threshold, where $H(x) = 1\ \forall \ x \ge x_0$ and $H(x) = 0\ \forall \ x<x_0$. The index $A_P$ was estimated through the retrievemodis function and stored in the object Parea (short for precipitation area):
## TRMM: daily data
retrievemodis() -> Parea
save(Parea,file='Parea.rda')
## data("Parea")
It was convenient to also include a number of other indices in Parea while estimating $A_P$, and the object also was designed to hold additional information such as the precipitation area over land and oceans respectively in addition to the precipiation intensity (wet-day mean precipitation) and total precipitation amount $P_{tot}=\int_A x da$ where $x$ is the precipitation. The total precipitation can also be written as $P_{tot}=A_P \overline{x}$.
attr(Parea,'longname')
#[1] "Precipitation area" "Precipitation area over ocean" "Precipitation area over land"
#[4] "wet-day mean preciptation" "total precipitation" "total precipitation over oceans"
#[7] "total precipitation over land" "number of valid data"
Quality check:
Repeat the calculations with NOAA/Pacific Marine Environment Lab (PMEL) Ferret:
```{bash, eval=FALSE}
use $1
set reg/y=50S:50
CANCEL LIST/HEAD
let precip = precipitation
let TP = precip[x=@din,y=@din]
let pre IF precip GT 1 then 1 else 0
let A = pre[x=@din,y=@din]
!save/append/clobber/file=totprecip.nc TP,A
list/append/file=totprecip.txt TP,A
cancel data ($1)
cancel memory
The bash-script `ferretjobs.sh` was used to carry out the ferret jobs:
```{bash, eval=FALSE}
## Ferret script to estimate the total precipitation and the area
## Rasmus.Benestad@physics.org
#!/bin/bash
for f in 3B42*.nc4; do
echo $f
/usr/local/ferret/bin/ferret -gif -script ~/TRMM/totprecip.jnl $f
done
The job is ran by typing the command:
```{bash,eval=FALSE}
./ferretjobs.sh
The results were then read in R and compared with the calculations done in R:
The data looks like
```{bash, eval=FALSE}
I / *: 104515. 8493.
I / *: 97743. 8408.
I / *: 98213. 8205.
I / *: 90036. 8106.
...
Which means that Ferret did not make use of the metadata (coordinate units etc) for properly estimating the area and the total precipitation.
tp.in <- read.table('~/TRMM/totprecip.txt')
pta <- cbind(tp.in$V4,tp.in$V5)
tpa.day <- zoo(pta,order.by=seq(from=as.Date('1998-01-01'),by='day',
length.out=length(pta[,1])))
names(tpa.day) <- c('tot.precip','precip.area')
data(Parea)
# Compare with the calculations done in R:
plot(stand(tpa.day[,2]),lwd=3,col='grey') # From PMEL Ferret
lines(stand(Parea[,1]),col='red',lty=3) # From R
## Total precipitation
plot(stand(tpa.day[,1]),lwd=3,col='grey')
lines(stand(Parea[,5]),col='red',lty=3)
The comparisons verified the consistency of the calulations in R and with Ferret, albeit with different units.
Repeating the analysis on monthly mean precipitation
A quality check was applied to the monthly TRMM data
## May have to run this several times with different values for yr1
## and having restarted R - there seems to be a memory leak in ncdf4.
TRMM2monthly()
## Processing for data("Parea.month")
The monthly TRMM files were also concatinated into one file using NCO so that the monthly data derived from daily $A_P$ could be compared with monthly data estimated from monthly mean precipitation:
```{bash, eval=FALSE}
ncrcat trmm-precip*.nc trmm-precip-mon.nc
```r
## Only the precipitation area - not total precipitation
monthP2area(x0=4.5) -> Parea.month
Quality check:
Repeat the calculations with PMEL Ferret:
```{bash, eval=FALSE}
use trmm-precip-mon.nc
set reg/y=50S:50
CANCEL LIST/HEAD
let TP = precip[x=@din,y=@din]
let pre IF precip GT 4.5 then 1 else 0
let A = pre[x=@din,y=@din]
save/clobber/file=totprecip.nc TP,A
list/clobber/file=totprecip.txt TP,A
The results were then read in R and compared with the calculations done in R:
```r
tp.in <- read.table('~/TRMM-monthly/totprecip.txt')
pta <- cbind(tp.in$V5,tp.in$V6)
tpa.mon <- zoo(pta,order.by=seq(from=as.Date('1998-01-01'),by='month',
length.out=length(pta[,1])))
names(tpa.mon) <- c('tot.precip','precip.area')
data(Parea.month)
# Compare with the calculations done in R:
plot(stand(tpa.mon[,2]),lwd=3,col='grey') # From PMEL Ferret
lines(stand(Parea.month[,1]),col='red',lty=3) # From R
The comparisons revealed some differences, as the two were not identical, possibly due to different ways of applying the threshold $x_0$ defining wet conditions.
ERAINT
ERAINT-based indices describing the global rate of evaporation $Q_s=\int_A Eda$, total atmospherc water, and the area with precipitation $A_P$.
## ERAINT
tp.in <- read.table('~/Downloads/totprecip.txt')
pta <- cbind(tp.in$V5,tp.in$V6,tp.in$V7,tp.in$V8,tp.in$V9)
n <- length(pta[,1])
PTA <- pta[seq(1,(n-1),by=2),] + pta[seq(2,n,by=2),]
tpa.eraint <- zoo(PTA,order.by=seq(from=as.Date('1979-01-01'),by='day',
length.out=length(PTA[,1])))
names(tpa.eraint) <- c('tot.precip','area.P.gt.0.1mm','area.P.gt.1mm','area.P.gt.10mm','tot.area')
save(tpa.eraint,file='tpa.eraint.rda')
ERAINT also suggest a decreasing $A_p$, from 36% to about 33%. This is roughly consistent with the TRMM data, but the total area estimated for ERAINT was estiamted to be $7.8 \times 10^8 km^2$ whereas the area of the earth is $5.1 \times 10^8 km^2$.
data(tpa.eraint)
plot(100*tpa.eraint[,2]/tpa.eraint[,5],col='grey',ylab=expression(A[P]))
grid()
lines(trend(100*tpa.eraint[,2]/tpa.eraint[,5]),lty=2)
Other quantities from ERAINT.
## Aggregated instantaneous surface moisture fluxes (evaporation)
moistflux() -> Qs
## data("Qs")
## Total atmospheric column water
columnH2O() -> H2O
## data("H2=")
From ERAINT, we also get the total columnt water: $H_2O=\int_A q da$ where $q$ is the water mass in the entire air column.
Alsothe total precipitation was estimated from ERAINT.
#ERAINT: units in m/s
prtot.eraint <- globalsum('~/Downloads/ERAINT_precip.nc')
## The data is stored as the daily mean precipitation for 00:00-12:00hr and 12:00-24:00hr bacthes
mm <- as.character(month(prtot.eraint))
mm[nchar(mm)==1] <- paste('0',mm[nchar(mm)==1],sep='')
yyyymm <- paste(year(prtot.eraint),mm,'01',sep='-')
prtot.eraint <- aggregate(prtot.eraint,yyyymm,mean)*1000 # Unsure about 'mean' or 'sum'
index(prtot.eraint) <- as.Date(index(prtot.eraint))
prtot.eraint <- prtot.eraint*1e3# units in mm rather than m
attr(prtot.eraint,'unit') <- rep('tons',3)
save(prtot.eraint,file='prtot.eraint.rda')
The units of ERAINT precipitation is $m$ and when aggregated using 'sum' the value is close to twice that of the TRMM data and twice the daily evaporation. The data comes with two values per month - it's not clear whether the both values are daily averages or 12-hour accumulated values.
Ferret script for extracting estimates of total precipitation totpre and precipitation area A bfrom ERAINT reanalysis between 50S-50N and on a monthle scale (using monthly precipitation):
```{bash, eval=FALSE}
ferret
yes? use ERAINT_precip.nc
yes? set reg/y=50S:50N
yes? CANCEL LIST/HEAD
yes? let totpre = TP[x=@din,y=@din]
yes? let pre IF TP GT 0.005 then 1 else 0
yes? let A = pre[x=@din,y=@din]
yes? list/clobber/file=totmonprecip.txt totpre,A
Ferret script for calculating total global evaporation from ERAINT reanalysis:
```{bash, eval=FALSE}
use ERAINT_Qs.nc
let Qs = IE[x=@din,y=@din]*86.4
CANCEL LIST/HEAD
list/file=qs.txt qs
R script to synthesise monthly data from two parts: 00:00-12:00 and 12:00-24:00 (this is the way the data was stored)
## TYhe ERAINT data were stored in a devious format:
#01-JAN-1979 00:00 / 1: 1.259E+12 1.504E+14
#01-JAN-1979 12:00 / 2: 1.278E+12 1.557E+14
#01-FEB-1979 00:00 / 3: 1.254E+12 1.573E+14
#01-FEB-1979 12:00 / 4: 1.275E+12 1.640E+14
tp.in <- read.table('~/Downloads/totmonprecip.txt')
pta <- cbind(tp.in$V5,tp.in$V6)
n <- length(pta[,1])
PTA <- pta[seq(1,(n-1),by=2),] + pta[seq(2,n,by=2),]
names(PTA) <- c('tot.precip','precip.area')
tpa.eraint <- zoo(PTA,order.by=seq(from=as.Date('1979-01-01'),by='month',
length.out=length(PTA[,1])))
Plotting the data (not shown)
data("prtot.eraint")
# Compare with the calculations done in R:
plot(stand(tpa.eraint[,1]),lwd=3,col='grey') # From PMEL Ferret
lines(stand(prtot.eraint[,1]),col='red',lty=3) # From R
print(mean(tpa.eraint[,1])*0.5)
print(mean(prtot.eraint[,1]))
The comparisons verified the consistency of the calulations in R and with Ferret, albeit with different units (not shown).
# Compare total precipitation between 50S-50N
plot(stand(tpa.mon[,1]),lwd=3,col='grey') # From TRMM
lines(stand(tpa.eraint[,1]),col='red',lty=3) # From ERAINT
Comparison for monthly precipitation area in ERAINT and TRMM (not shown)
# Compare precipitation area from monthly mean precipitation
plot(stand(tpa.mon[,2]),lwd=3,col='grey') # From TRMM
lines(stand(tpa.eraint[,2]),col='red',lty=3) # From ERAINT
CMIP 5 simulations
Estimate monthly aggregated precipitation area from the RCP4.5 CMIP5 experiment. These are only crude indices since the spatial resolution varies from model to model and the grid boxes represent an area average of the rainfall while the rain in reality only falls over a fraction of the area.
PareaCMIP() -> Parea.cmip5
Total precipitation
demo('global.t2m.gcm',ask=FALSE)
ptot.cmip <- globalmean(path='~/CMIP5.monthly/rcp45',
annual=TRUE,pattern='pr_',param='pr',anomaly=FALSE,FUN='sum')
save(ptot.cmip,file='ptot.cmip.rda')
Loading the data for the analysis
For the analysis, the data/indices need to be loaded in memory:
data(Parea)
data(Parea.month) # Only precipitation area
data(Parea.eraint.month) # Only precipitation area
data(prtot.eraint)
Some basic statistics such as the mean area of precipitation (km$^2$), wet-day mean precipitation (mm/day) and total precipitation (Ktons) estimated for the 50S-50N according to the TRMM data:
colMeans(Parea)
The analysis
Figures 1-4.
This code is provided for making the analysis more transparent and easier to replicate. However, the plotsa are shown in the main paper.
Fig 1 shows a map of a random day of precipirtation just to illustrate what is meant by the precipitation area (blue regions).
Fig1()
Fig 2 shows the daily precipitation area estimated from TRMM data.
Fig2()
t1 <- start(Parea)
t2 <- end(Parea)
A.p <- window(trend(Parea[,1]),start=t1,end=t2)
## Lower-Upper range
print(range(Parea[,1]))
## P-value for trend analysis
print(trend.pval(Parea[,1]))
## Trend: min and max
print(range(A.p))
## Change in percentage
print(100*diff(range(A.p))/max(A.p))
## Trend as percentage of the total area:
print(round(100*range(A.p)/attr(Parea,'total area'),2))
Fig 3 compares monthly mean area estimated from daily area (TRMM and ERAINT) and monthly area estimated from monthly precipitation TRMM and ERAINT). Also compare monthly area between 50S--50N with global value. Fig 3 shows the daily precipitation area indices: both the daily indices and $A_P$ estimated from monthly mean precipitation for both. TRMM and ERAINT. The monthly aggregate differs from the daily area which seems to be shifted geographically from day to day, e.g. through storm tracks.
Fig3()
The monthly mean of daily indices do not follow the indeces of monthly mean precipitation. It is not expected that they should correspond, as they represent different aspects.
The time scale of preciptation is approximately hours, and often the precipitation does not move very far over the course of a day. Nevertheless, the 24-hr precipitation also involves moving storms and weather fronts to some extent, and are likely to over-estimate the area of precipitation.
The area of monthly mean precipitation is expected to be affected by the advection and migration of phenomena producing precipitation. For instance, the daily precipitation area $A_P$ could be constant, but the precipitation phenomena may move around so that the area of the monthly mean precipitation would appear to be larger. Even when setting a threshold $x_0$ for defining precipitation, the longer time scale will blur the area.
The spatial resolution too has an effect on the estimated area for a similar reason as the longer time scales. It is expected that it only rains in a fraction of the grid box, and the higher spatial resolution, the more precise we can expect the results to be.
The threshold $x_0$ was set to 1mm/day for calculation of the daily $A_P$but $x_0 = 4.5 mm/day$ for $A_P$ derived from monthly mean precipitation, except for ERANIT which had a coarser spatial resolution.
Fig 4 presents precipitation area indices for the CMIP5 RCP4.5 simulations were estimated from monthly mean precipitation $A_P = \sum_i a_i P_i\ \forall P_i > x_0$. The colour of the curves indicate their spatial resolution marked by the legend (number represents degree in logitude; red higher and blue lower spaital resolution).
Fig4()
Fig 4 also provides a summary of the trends indicated by the GCMs (%/decade). The mean increase over 2000--2100 is 2.4%.
Precipitation area indices for the CMIP5 RCP4.5 simulations were based on monthly mean precipitation. The CMIP5 results suggest an increase in the precipitation area for the monthly mean precipitation. Also the TRMM-based precipitation area based on monthly mean precipitation hinted towards an increase over time, suggesting more variable rainfall patterns given that the daily $A_P$ has declined. However, it is uncertain whether the GCMs also suggest a reduced $A_P$ for daily precipitation in the future.
For the CMIP5 data, $x_0 = 4 mm/day$ for calculating $A_P$.
Table GCM and mean precipitation area
Check how the mean precipitation area in (1000km)^2 varies wih the GCM and spatial resolution:
data(Parea.cmip5)
y <- colMeans(Parea.cmip5)*1.e-6
names(y) <- attr(Parea.cmip5,'model_id')
print(y)
Test to see if there is a systematic dependency on the model spatial resolution.
## ANOVA to analyse potential systematic relationship between model resolution and mean precipitation area
data("IPCC.AR5.Table.9.A.1")
nm1 <- tolower(gsub('.','',gsub('-','',IPCC.AR5.Table.9.A.1$Model.Name),fixed=TRUE))
nm2 <- tolower(gsub('.','',gsub('-','',attr(Parea.cmip5,'model_id')),fixed=TRUE))
##x <- IPCC.AR5.Table.9.A.1$Horizontal.Grid[charmatch(nm2,nm1)]
x <- cmipgcmresolution()[charmatch(nm2,nm1)]
rbind(IPCC.AR5.Table.9.A.1$Model.Name,as.character(IPCC.AR5.Table.9.A.1$Horizontal.Grid),x)
resbias <- data.frame(y=y,x=x)
print(summary(lm(y ~ x, data=resbias)))
The results from the ANOVA suggests that there is a tendency for smaller mean precipitation area for models with higher resolution since the regression coefficient $x$ is negative. The relationship is statistically significant at the 5\%-level. One issue is that different spatial resolution may mean different values, as they reflect area mean values for the grid box. This may also influence the analysis as a constant threshold value $x_0$ for all these models may select precipitation events differently for models with high and low spatial resolution, unless the regridding of the models has taken the area-mean into account.
The resolution in this case was derived from a table of model details from the IPCC AR5, and the data here had all been regridded to the same spatial grid in the KNMI Climate Explorer https://climexp.knmi.nl/.
Supporting material and additional analysis
Rainfall area over ocean
The precipitation area over ocean was esimated by masking out all land/contentent grid boxes. This index is the third in the Parea object:
timeseries(Parea,is=2)
print(trend.coef(Parea[,2]))
print(trend.pval(Parea[,2]))
print(range(window(trend(Parea[,2]),start=t1,end=t2)))
There was a clear decrease in $A_P$ over the oceans.
Rainfall area over land
Consequently, the rainfall area over land was estimated by masking all ocean grid points:
timeseries(Parea,is=3)
print(trend.coef(Parea[,3]))
print(trend.pval(Parea[,3]))
print(range(window(trend(Parea[,3]),start=t1,end=t2)))
The total rain area over land is minor compared to the rain area over the oceans. Furthermore, the trend is much weaker and not statistically significant at the 5\%-level. There is a more pronounced annual cycle in $A_P$ over land than over the oceans, which may be related to monsoon-type systems.
Total precipitation from the TRMM data
The 50S-50N total precipitation amount was estimated by summing up all grid boxes with rainfall for each day:
timeseries(Parea,is=5,denominator=1e6,ylab='Gton')
## Lower-Upper range
print(range(Parea[,5])/1.0e6)
## Check period
print(c(start(Parea),end(Parea)))
## P-value for trend analysis
print(trend.coef(Parea[,5]))
print(trend.pval(Parea[,5]))
## Trend: min and max: ocean and land
print(range(window(trend(Parea[,5]),start=t1,end=t2)/1e6))
## Trend: min and max: ocean
print(trend.coef(Parea[,6]))
print(trend.pval(Parea[,6]))
print(range(window(trend(Parea[,6]),start=t1,end=t2)/1e6))
## Trend: min and max: land
print(trend.coef(Parea[,7]))
print(trend.pval(Parea[,7]))
print(range(window(trend(Parea[,7]),start=t1,end=t2)/1e6))
The total daily volume of rainfall between 50S--50N exhibits a wide range of daily fluctuations between 868 -- 1494 Gton, but there is no strong long-term trend. An increasing trend is seen in the plot from 1120 -- 1150 Gton.
The total precipitation from ERAINT:
## ERAINT
plot(prtot.eraint[,1]/1e9,ylab='Gton',main='ERAINT total precipitation 50S-50N',col='red',lwd=2)
lines(trend(prtot.eraint[,1]/1e9))
## Trend: min and max
print(trend.coef(prtot.eraint[,1]))
print(trend.pval(prtot.eraint[,1]))
print(range(window(trend(prtot.eraint[,1]),start=t1,end=t2)/1e9))
Compare the total precipitation between 50S--50N with the rest of the globe:
print(colMeans(prtot.eraint)/1.0e9)
Rainfall amount over ocean
This index can be split into the total rainfall amount over the oceans by masking out the land grid boxes (not shown)
## The units are mm*km^2 -> 1000m^3 -> million kg or 1000 tons (1 Kton).
timeseries(Parea,is=6,denominator=1e6,ylab='Gton')
The contribution to the total precipitation from the oceans exhibits a slightly more pronounced increasing trend than the total oceans+land amounts do.
The fraction of earth's area between 50S--50N is $\int$ $2 \pi a \cos (\phi) a d \phi/4 \pi a^2 = [\sin (\phi_2) - \sin (\phi_1)]/2$ where $\phi= \pm \pi*50/180$.
print(paste(round(100*sin(pi*50/180)),"% of earth's surface area is between 50S-50N"))
Rainfall amount over land
The corresponding index for the total land precipitation is (not shown)
timeseries(Parea,is=7,denominator=1.0e6,ylab='Gton')
The TRMM data suggests a declining trend in the total precipitation amount land in contrast to the precipitation over ocean and the sum oceans+land.
Precipitation intensity from TRMM
The rain intensity $\mu$ can be estimated by taking the average of all grid boxes with rain for each day:
timeseries(Parea,is=4,denominator=1)
## Trend: min and max
print(trend.coef(Parea[,4]))
print(trend.pval(Parea[,4]))
print(round(range(window(trend(Parea[,4]),start=t1,end=t2)),2))
The graph reveals an increasing trend in $\overline{P}$, and the trend analysis is applied to a long series of daily values which makes it statistical significant at the 5\%-level (p-value $\approx 10^{-21}$).
A combination of a declining rainfall area over time with increasing precipitation amount is consistent with an increasing precipitation intensity: $P = \overline{P} * A_P$.
Check - covariation with the global mean temperature?
Assess the scaling relation between observed precipitation area and the global mean temperature through regression analysis ($A_P = \beta_0 + \beta_1 T + \eta$) where regression coefficient $\beta_1$ describes the slope (rate of change in $A_P$ with the global mean temperature) and $\beta_0$ the zero-intersect and $\eta$ being a nooise-term.
require(esd)
hc4 <- HadCRUT4()
## The precipitation area:
XY <- merge(annual(hc4),annual(Parea[,1]),all=FALSE)
cal <- data.frame(x=coredata(XY[,1]),y=coredata(XY[,2]))
pax <- lm(y ~x, data=cal)
print(summary(pax))
## Total precipitation:
XY <- merge(annual(hc4),annual(Parea[,5]),all=FALSE)
cal <- data.frame(x=coredata(XY[,1]),y=coredata(XY[,2]))
print(summary(lm(y ~x, data=cal)))
An ANOVA on the TRMM annual mean values indicates strong covariance between the global mean temperature and the precipitation area, with $R^2$ of 0.6 and a negative regression coefficient $\beta_1 = -18 \times 10^6 km^2/$degree K. (statistically sign. at the 0.1\%-level). A similar analysis for the total precipitation suggests little covariance (not shown).
## Compare annual cycles
xy <- merge(stand(window(hc4,start=t1,end=t2)),stand(Parea[,1]),all=TRUE)
plot(aggregate(xy,month,mean,na.rm=TRUE),col=c('blue','red'),plot.type='single')
The there is little covariance between the mean seasonal cycles of HadCRUT4 and $P_A$.
Naive extrapolation for future daily $A_P$ based on global mean temperature change
data("global.t2m.cmip5")
Pacmip5 <- global.t2m.cmip5
for (i in 1:dim(global.t2m.cmip5)[2]) {
pre <- data.frame(x=coredata(global.t2m.cmip5[,i]))
coredata(Pacmip5)[,i] <- predict(pax,newdata=pre)
}
plot(Pacmip5*1.0e-6,plot.type='single',main='Extrapolated precipitation area',
ylab='million km^2',ylim=c(0,200),col=rgb(0,0,0,0.2))
grid()
print(100*(mean(window(Pacmip5,start=2079,end=2099)) - mean(window(Pacmip5,start=1998,end=2016))) /mean(window(Pacmip5,start=1998,end=2016)))
A naive extrapolation assuming that (a) the regression analysis between the global mean temperature and precipitation area captures a true dependency and (b) that this dependency is stationary and holds in the future suggests a reduction in $A_P$ from near 110 million km$^2$ in the 20th century to 75 million km$^2$ by 2100. The CMIP5 ensemble mean suggests around a 28\% decrease in $A_P$.
Total atmospheric water content from ERAINT
The total 50S-50N atmospheric water content was estimated from the ERAINT reanalysis
data(H2O)
plot(H2O[,1]/1e9,ylab='Gton')
## The 50S-50N water content is much higher than for the high latitude regions
## "kg km**-2" agregated (sum) over the area
lines(trend(H2O[,1]/1e9))
grid()
## Trend: min and max
print(trend.coef(H2O[,1]))
print(trend.pval(H2O[,1]))
print(range(window(trend(H2O[,1]),start=t1,end=t2)/1e9))
print(c(start(H2O),end(H2O)))
colMeans(H2O)/1e9 ## Gtons
The total atmospheric water exhibits weak long-term changes according to the ERAINT reanalysis and shows a hint of an increase over time.
The surface moisture flux $Q_s$ from ERAINT is the same as evaporation $E_s$ and the declining trends (negative numbers) indicates an increased evaporation over time. This is consistent with increased $P_{tot}$ in the TRMM data. Negative values indicate flux out of the oceans.
data(Qs)
plot(annual(Qs,FUN='mean'),ylab='Gton')
## The 50S-50N moisture flux is much higher than for the high latitude regions
## Units are mm/day * km^2 which is equivalent to 1000m^3/day or kilo-tons/day
lines(trend(annual(Qs,FUN='mean')))
grid()
## Giga-tons/day - divide by 1e6
## Trend: min and max
print(trend.coef(Qs[,1]))
print(trend.pval(Qs[,1]))
print(range(window(trend(Qs[,1]),start=t1,end=t2)))
print(c(start(Qs),end(Qs)))
mean(Qs)
Table on the amount of evaporation and precipitation - water fluxes
Table 1:
print('Moisture flux (ERAINT; Gton/day)')
print(round(mean(Qs)))
print('Column water (ERAINT; Gton)')
print(round(colMeans(H2O)*1.e-9))
print('Precipitation (ERAINT; Gton)')
print(round(colMeans(prtot.eraint)*1.e-9))
print('Mean total precip area 50S-50N (Gton/day)')
print(round(mean(Parea[,1]*1.e-6)))
print('Mean total precip area over ocean 50S-50N (Gton/day)')
print(round(mean(Parea[,2]*1.e-6)))
print('Mean total precip area 50S-50N over land (Gton/day)')
print(round(mean(Parea[,3]*1.e-6)))
print('Mean total precip amount 50S-50N (Gton/day)')
print(round(mean(Parea[,5]*1.e-6)))
print('Mean total precip amount over ocean 50S-50N (Gton/day)')
print(round(mean(Parea[,6]*1.e-6)))
print('Mean total precip amount 50S-50N over land (Gton/day)')
print(round(mean(Parea[,7]*1.e-6)))
Table 2:
h2o <- window(trend(H2O[,1]),start=t1,end=t2)/1e15
print(round(c(h2o[1],h2o[length(h2o)])))
qs <- window(trend(Qs),start=t1,end=t2)
print(round(c(qs[1],qs[length(qs)])))
y <- window(trend(Parea[,1]),start=t1,end=t2)/1e6
print(round(c(y[1],y[length(y)])))
y <- 100*window(trend(Parea[,1]),start=t1,end=t2)/attr(Parea,'total area')
print(round(c(y[1],y[length(y)])))
y <- window(trend(Parea[,2]),start=t1,end=t2)/1e6
print(round(c(y[1],y[length(y)])))
y <- window(trend(Parea[,3]),start=t1,end=t2)/1e6
print(round(c(y[1],y[length(y)])))
y <- window(trend(Parea[,5]),start=t1,end=t2)/1e6
print(round(c(y[1],y[length(y)])))
y <- window(trend(prtot.eraint[,1]),start=t1,end=t2)/1e9
print(round(c(y[1],y[length(y)])))
y <- window(trend(Parea[,6]),start=t1,end=t2)/1e6
print(round(c(y[1],y[length(y)])))
y <- window(trend(Parea[,7]),start=t1,end=t2)/1e6
print(round(c(y[1],y[length(y)])))
y <- window(trend(Parea[,4]),start=t1,end=t2)
print(round(c(y[1],y[length(y)]),1))
y <- window(trend(Parea.month[,1])/1e6,start=t1,end=t2)
print(round(c(y[1],y[length(y)]),1))
y <- window(trend(tpa.eraint[,3]),start=t1,end=t2)/1e12 # Average of two
print(round(c(y[1],y[length(y)]),1))
index(Parea.eraint.month) <- as.Date(index(Parea.eraint.month))
y <- window(trend(Parea.eraint.month),start=t1,end=t2)/1e6
print(round(c(y[1],y[length(y)]),1))
Previous estimates of global precipitation flux
Table 1 in Zektser et al (1993) http://www.sciencedirect.com/science/article/pii/0022169493901829
## Example in the text: evaporation
print(294*128) ## km^3/year - should be 38,000 km^3 according to the article
## As kg/day: 1 km^3 = 1.oe9 tons = 1 Gton.
print(294*128/365.25)
## Precipitation over land
print(paste('Precipitation: 834 mm/year=',round(834*128/365.25),'Gton/day'))
Estimate of the global mean precipitation from Kiel and Trenberth: 2.46 -- 2.90 mm/day http://journals.ametsoc.org/doi/pdf/10.1175/1520-0477%281997%29078%3C0197%3AEAGMEB%3E2.0.CO%3B2
print(paste('2.46mm/day =',round(2.46*1.0e-3 * 4*pi*6378000^2*1.0e-9),'Gton/day'))
print(paste('2.46mm/day =',round(2.90*1.0e-3 * 4*pi*6378000^2*1.0e-9),'Gton/day'))
Legates (1995): 2.6 to 3.1 mm/day http://onlinelibrary.wiley.com/doi/10.1002/joc.3370150302/full
print(paste('2.6mm/day =',round(2.6*1.0e-3 * 4*pi*6378000^2*1.0e-9),'Gton/day'))
print(paste('3.1mm/day =',round(3.1*1.0e-3 * 4*pi*6378000^2*1.0e-9),'Gton/day'))
Test the correlation between the precipitation area and galactic cosmic rays
Compare the cosmic rays from Climax with the precipitation area (not shown)
require(replicationDemos)
data(gcr)
t <- paste(substr(gcr$Date,7,10),substr(gcr$Date,1,2),substr(gcr$Date,4,5),sep='-')
climax <- zoo(as.numeric(as.character(gcr$Climax)),order.by=as.Date(t))
plot(climax,Parea[,1])
gcrpa <- merge(climax,Parea[,1],all=FALSE)
ok <- is.finite(gcrpa[,1]) & is.finite(gcrpa[,2])
print(cor.test(gcrpa[ok,1],gcrpa[ok,2]))
Sanity check - test of consistency
Compare the total precipitation from the daily data and the monthly data
P.tot <- merge(aggregate(Parea[,5],year,mean),
aggregate(prtot.eraint[,1]/1000,year,mean),all=TRUE)
plot(P.tot/1.0e6,plot.type='single',col=c('black','red'),led=c(2,2))
Estimate of the portion of the vertical energy flow is due to convection from the total precipitation
The precipitation was assumed to be a by-product of moist convection with an energy content equivalent to the release of latent heat of evaporation $P L_E$:
$f_{conv}= \frac{P L_e}{\pi a^2 (1-A) S_0}$
a <- 6378000 # Earth's area (m)
A = 0.3 # The albedo
S0 <- 1361 # The Solar constant (W/m²)
Le <- 2264.76 * 1000 # Latent heat of evaporation (https://en.wikipedia.org/wiki/Latent_heat) J/kg
tot.E <- pi * a^2 * (1-A) * S0 * 24 * 3600 # Energy over 24 hrs
P <- sum(colMeans(prtot.eraint))*1000
frac.E <- P*Le/tot.E * 100
print(paste('Convection is responsible for ',round(frac.E,2),'percent of the vertical heat flow'))
Trend anlysis of the proportn of vertical energy flow being connected to convection
fE <- (prtot.eraint[,1]+prtot.eraint[,2]+prtot.eraint[,3])*Le*100000/tot.E
plot(fE)
grid()
lines(trend(fE))
prop <- window(trend(fE),start=t1,end=t2)
print(round(c(prop[1],prop[length(prop)]),1))
For the TRMM data:
y <- window(trend(Parea[,5]),start=t1,end=t2)
print(paste('Percentage energy flow associated with convection',round(c(y[1],y[length(y)])*Le*1e6/tot.E*100),'%'))
Estimate the globally summed precipitation from CMIP5:
data(ptot.cmip)
plot(zoo(ptot.cmip)/1e6,plot.type='single',ylab='Gton')
pclim <- colMeans(window(ptot.cmip,start=1980,end=2009))/1e6 # mean total precipitation
print(summary(pclim*Le/tot.E * 100))
print(summary(pclim)) # In gigatons
Ptrends <- apply(window(ptot.cmip,start=2000,end=2100)/1e6,2,trend.coef)
print(summary(Ptrends))
Ftrends <- apply(100*window(ptot.cmip,start=2000,end=2100)*Le*1e6/tot.E,2,trend.coef)
## Change in percent of energy flow connected with convection between 2000 and 2100
print(summary(Ftrends*10))
Acknowledgements:
The PMEL/NOAA software ferret http://ferret.pmel.noaa.gov/Ferret/ and David Pierce's ncview were used to check and verify the netCDF files and the calculations. The function globalsumand ferret indicated similar results for the monthly mean precipitation totals.
brasmus/preciparea documentation built on May 7, 2019, 11:11 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
Statistics on the precipitation area
The instantaneous area of precipitation is an important indicator describing the hydrological cycle and the state of earth's climate. This is a document of data analysis carried out to elucidate this aspect of the climate, based on model data and historical observations, and is intended as the supporting material for the paper with the title 'Was the 1998--2016 decrease in the tropical precipitation area real?'. It is meant to provide a record like a lab note-book, to enhance the possibility of replication and provide details about the analysis. It provides openness and transparency, and is not mean to be a 'well-structured' document.
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
Getting the data
The TRMM data was downloaded (February 2017) from NASA data portals https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/ using he instructions provided in https://disc.sci.gsfc.nasa.gov/recipes/?q=recipes/How-to-Download-Data-Files-from-HTTP-Service-with-wget and https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/doc/TRMM_Readme_v3.pdf. It is necessary to create a user account and set up user details/coockies as descibed in the intructions. The code presented here was written for Linux, and it is uncertain whether it works for Windows machines. It maybe will work on iOS. The data volumes are large and require some time and space for a complete download
Satellite data - the Tropical Rain Measurement Mission (TRMM)
Download the daily TRMM data (~16Gb and takes time) on a Linux platform (using a terminal and bash-scripts).
```{bash, eval=FALSE}
Retrieving the data:
https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/
See this web page for how to download the files.
https://disc.sci.gsfc.nasa.gov/recipes/?q=recipes/How-to-Download-Data-Files-from-HTTP-Service-with-wget
https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/doc/TRMM_Readme_v3.pdf
Script for downloading thr TRMM data
!/bin/bash
mkdir TRMM cd TRMM for year in {1998..2016} do for month in {1..12} do wget --load-cookies ~/.urs_cookies --save-cookies ~/.urs_cookies --keep-session-cookies -r -c -nH -nd -np -A nc4 "https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_L3/TRMM_3B42_Daily.7/$year/$month" done done
#### Independent reanalysis - ERAINT
The ERAINT data were downloaded from <http://apps.ecmwf.int/datasets/data/interim-full-moda/levtype=sfc/> (Monthly means of Daily means and daily precipitation) as netCDF and renamed `ERAINT_t2m.nc` (2-meter temperature), `ERAINT_Qs.nc` (surface moisture flux), and `ERAINT_Q-columntotal.nc` (total column water). Te daily precipitation was downloaded in multiple files due to large data volumes: `ERAINT_24hr_precip*.nc`. The unit of the precipitation was 'm' (it is assumed it was m/day).
The following lines of code were the contents of the script `totprecip.jnl` implemnted by PMEL Ferret:
```{bash, eval=FALSE}
use $1
set reg/y=50S:50
CANCEL LIST/HEAD
let TotP = TP[x=@din,y=@din]
let pre001 IF TP GT 0.0001 then 1 else 0
let Ap001 = pre001[x=@din,y=@din]
let pre01 IF TP GT 0.001 then 1 else 0
let Ap01 = pre01[x=@din,y=@din]
let pre1 IF TP GT 0.01 then 1 else 0
let Ap1 = pre1[x=@din,y=@din]
let pre0 IF TP LT -9999 then 0 else 1
let TotA=pre0[x=@din,y=@din]
!save/append/clobber/file=totprecip.nc TotP,Ap01,Ap1,TotA
list/append/file=totprecip.txt TotP, Ap001,Ap01,Ap1,TotA
cancel data ($1)
cancel memory
! Earth's area = 5.111859e+14 m^2
The scripty was ran through a bash job:
```{bash, eval=FALSE}
!/bin/bash
for f in ERAINT_24hr-precip*.nc; do echo $f /home/rasmusb/ferret/bin/ferret -gif -script ~/Downloads/totprecip.jnl $f done
The use of MERRA-2 was also attempted, but without success. The right data (daily precipitation) was hard to find, and a search on the web site gave numerous hits with all sorts of precipitation results but not for 24-hr precipitation on a daily basis through the reanalysis interval. A Google search did not give any useful tips. An email request was never answered. ### The R-code To keep all R-scripts/functions, extracted data and documentation together, a dedicated R-package `preciparea` was created that must be installed and activated for replicating this work. It is available at <https://github.com/brasmus/preciparea>.The R-scripts are open-code in this package and is part of the supporting material (SM) and contains both R-scripts and extracted data for this analysis. The R-package can be installed from github: ```r install.esd <- ("esd" %in% rownames(installed.packages()) == FALSE) install.preciparea <- ("preciparea" %in% rownames(installed.packages()) == FALSE) install.devtools <- ("devtools" %in% rownames(installed.packages()) == FALSE) & (install.esd | install.preciparea) if (install.devtools) { print('Need to install the devtools package') ## You need online access. install.packages('devtools') } library(devtools) if (install.esd) install_github('metno/esd') if (install.preciparea) install_github('brasmus/preciparea')
The package must be activated in R:
library(preciparea) library(esd) sessionInfo()$otherPkgs
Preparing indices from the data
This section shows how indeces stored in the preciparea package were derived, such as data("Parea").
TRMM
The TRMM netCDF files were unfortunately not stored with CF convention and a time dimension, and it was unclear how common netCDf processing tools like ncra could be applied to estimate the monthly means. Instead, the monthly data were aggrigated from the daily ones in the R-environment and the R-package preciparea (which was time-consuming).
The following command line were used to prepare the daily indicex describing the precipitation area $A_P = \int_A H(x-x_0)da$, where $H$ is the Heaviside function, $x$ is the precipitation, and $x_0$ is a precipitation threshold, where $H(x) = 1\ \forall \ x \ge x_0$ and $H(x) = 0\ \forall \ x<x_0$. The index $A_P$ was estimated through the retrievemodis function and stored in the object Parea (short for precipitation area):
## TRMM: daily data retrievemodis() -> Parea save(Parea,file='Parea.rda') ## data("Parea")
It was convenient to also include a number of other indices in Parea while estimating $A_P$, and the object also was designed to hold additional information such as the precipitation area over land and oceans respectively in addition to the precipiation intensity (wet-day mean precipitation) and total precipitation amount $P_{tot}=\int_A x da$ where $x$ is the precipitation. The total precipitation can also be written as $P_{tot}=A_P \overline{x}$.
attr(Parea,'longname') #[1] "Precipitation area" "Precipitation area over ocean" "Precipitation area over land" #[4] "wet-day mean preciptation" "total precipitation" "total precipitation over oceans" #[7] "total precipitation over land" "number of valid data"
Quality check:
Repeat the calculations with NOAA/Pacific Marine Environment Lab (PMEL) Ferret:
```{bash, eval=FALSE}
use $1
set reg/y=50S:50
CANCEL LIST/HEAD
let precip = precipitation
let TP = precip[x=@din,y=@din]
let pre IF precip GT 1 then 1 else 0
let A = pre[x=@din,y=@din]
!save/append/clobber/file=totprecip.nc TP,A
list/append/file=totprecip.txt TP,A
cancel data ($1)
cancel memory
The bash-script `ferretjobs.sh` was used to carry out the ferret jobs:
```{bash, eval=FALSE}
## Ferret script to estimate the total precipitation and the area
## Rasmus.Benestad@physics.org
#!/bin/bash
for f in 3B42*.nc4; do
echo $f
/usr/local/ferret/bin/ferret -gif -script ~/TRMM/totprecip.jnl $f
done
The job is ran by typing the command:
```{bash,eval=FALSE} ./ferretjobs.sh
The results were then read in R and compared with the calculations done in R:
The data looks like
```{bash, eval=FALSE}
I / *: 104515. 8493.
I / *: 97743. 8408.
I / *: 98213. 8205.
I / *: 90036. 8106.
...
Which means that Ferret did not make use of the metadata (coordinate units etc) for properly estimating the area and the total precipitation.
tp.in <- read.table('~/TRMM/totprecip.txt') pta <- cbind(tp.in$V4,tp.in$V5) tpa.day <- zoo(pta,order.by=seq(from=as.Date('1998-01-01'),by='day', length.out=length(pta[,1]))) names(tpa.day) <- c('tot.precip','precip.area')
data(Parea) # Compare with the calculations done in R: plot(stand(tpa.day[,2]),lwd=3,col='grey') # From PMEL Ferret lines(stand(Parea[,1]),col='red',lty=3) # From R
## Total precipitation plot(stand(tpa.day[,1]),lwd=3,col='grey') lines(stand(Parea[,5]),col='red',lty=3)
The comparisons verified the consistency of the calulations in R and with Ferret, albeit with different units.
Repeating the analysis on monthly mean precipitation
A quality check was applied to the monthly TRMM data
## May have to run this several times with different values for yr1 ## and having restarted R - there seems to be a memory leak in ncdf4. TRMM2monthly() ## Processing for data("Parea.month")
The monthly TRMM files were also concatinated into one file using NCO so that the monthly data derived from daily $A_P$ could be compared with monthly data estimated from monthly mean precipitation:
```{bash, eval=FALSE} ncrcat trmm-precip*.nc trmm-precip-mon.nc
```r ## Only the precipitation area - not total precipitation monthP2area(x0=4.5) -> Parea.month
Quality check:
Repeat the calculations with PMEL Ferret:
```{bash, eval=FALSE} use trmm-precip-mon.nc set reg/y=50S:50 CANCEL LIST/HEAD let TP = precip[x=@din,y=@din] let pre IF precip GT 4.5 then 1 else 0 let A = pre[x=@din,y=@din] save/clobber/file=totprecip.nc TP,A list/clobber/file=totprecip.txt TP,A
The results were then read in R and compared with the calculations done in R: ```r tp.in <- read.table('~/TRMM-monthly/totprecip.txt') pta <- cbind(tp.in$V5,tp.in$V6) tpa.mon <- zoo(pta,order.by=seq(from=as.Date('1998-01-01'),by='month', length.out=length(pta[,1]))) names(tpa.mon) <- c('tot.precip','precip.area')
data(Parea.month) # Compare with the calculations done in R: plot(stand(tpa.mon[,2]),lwd=3,col='grey') # From PMEL Ferret lines(stand(Parea.month[,1]),col='red',lty=3) # From R
The comparisons revealed some differences, as the two were not identical, possibly due to different ways of applying the threshold $x_0$ defining wet conditions.
ERAINT
ERAINT-based indices describing the global rate of evaporation $Q_s=\int_A Eda$, total atmospherc water, and the area with precipitation $A_P$.
## ERAINT tp.in <- read.table('~/Downloads/totprecip.txt') pta <- cbind(tp.in$V5,tp.in$V6,tp.in$V7,tp.in$V8,tp.in$V9) n <- length(pta[,1]) PTA <- pta[seq(1,(n-1),by=2),] + pta[seq(2,n,by=2),] tpa.eraint <- zoo(PTA,order.by=seq(from=as.Date('1979-01-01'),by='day', length.out=length(PTA[,1]))) names(tpa.eraint) <- c('tot.precip','area.P.gt.0.1mm','area.P.gt.1mm','area.P.gt.10mm','tot.area') save(tpa.eraint,file='tpa.eraint.rda')
ERAINT also suggest a decreasing $A_p$, from 36% to about 33%. This is roughly consistent with the TRMM data, but the total area estimated for ERAINT was estiamted to be $7.8 \times 10^8 km^2$ whereas the area of the earth is $5.1 \times 10^8 km^2$.
data(tpa.eraint) plot(100*tpa.eraint[,2]/tpa.eraint[,5],col='grey',ylab=expression(A[P])) grid() lines(trend(100*tpa.eraint[,2]/tpa.eraint[,5]),lty=2)
Other quantities from ERAINT.
## Aggregated instantaneous surface moisture fluxes (evaporation) moistflux() -> Qs ## data("Qs") ## Total atmospheric column water columnH2O() -> H2O ## data("H2=")
From ERAINT, we also get the total columnt water: $H_2O=\int_A q da$ where $q$ is the water mass in the entire air column.
Alsothe total precipitation was estimated from ERAINT.
#ERAINT: units in m/s prtot.eraint <- globalsum('~/Downloads/ERAINT_precip.nc') ## The data is stored as the daily mean precipitation for 00:00-12:00hr and 12:00-24:00hr bacthes mm <- as.character(month(prtot.eraint)) mm[nchar(mm)==1] <- paste('0',mm[nchar(mm)==1],sep='') yyyymm <- paste(year(prtot.eraint),mm,'01',sep='-') prtot.eraint <- aggregate(prtot.eraint,yyyymm,mean)*1000 # Unsure about 'mean' or 'sum' index(prtot.eraint) <- as.Date(index(prtot.eraint)) prtot.eraint <- prtot.eraint*1e3# units in mm rather than m attr(prtot.eraint,'unit') <- rep('tons',3) save(prtot.eraint,file='prtot.eraint.rda')
The units of ERAINT precipitation is $m$ and when aggregated using 'sum' the value is close to twice that of the TRMM data and twice the daily evaporation. The data comes with two values per month - it's not clear whether the both values are daily averages or 12-hour accumulated values.
Ferret script for extracting estimates of total precipitation totpre and precipitation area A bfrom ERAINT reanalysis between 50S-50N and on a monthle scale (using monthly precipitation):
```{bash, eval=FALSE} ferret yes? use ERAINT_precip.nc yes? set reg/y=50S:50N yes? CANCEL LIST/HEAD yes? let totpre = TP[x=@din,y=@din] yes? let pre IF TP GT 0.005 then 1 else 0 yes? let A = pre[x=@din,y=@din] yes? list/clobber/file=totmonprecip.txt totpre,A
Ferret script for calculating total global evaporation from ERAINT reanalysis:
```{bash, eval=FALSE}
use ERAINT_Qs.nc
let Qs = IE[x=@din,y=@din]*86.4
CANCEL LIST/HEAD
list/file=qs.txt qs
R script to synthesise monthly data from two parts: 00:00-12:00 and 12:00-24:00 (this is the way the data was stored)
## TYhe ERAINT data were stored in a devious format: #01-JAN-1979 00:00 / 1: 1.259E+12 1.504E+14 #01-JAN-1979 12:00 / 2: 1.278E+12 1.557E+14 #01-FEB-1979 00:00 / 3: 1.254E+12 1.573E+14 #01-FEB-1979 12:00 / 4: 1.275E+12 1.640E+14 tp.in <- read.table('~/Downloads/totmonprecip.txt') pta <- cbind(tp.in$V5,tp.in$V6) n <- length(pta[,1]) PTA <- pta[seq(1,(n-1),by=2),] + pta[seq(2,n,by=2),] names(PTA) <- c('tot.precip','precip.area') tpa.eraint <- zoo(PTA,order.by=seq(from=as.Date('1979-01-01'),by='month', length.out=length(PTA[,1])))
Plotting the data (not shown)
data("prtot.eraint") # Compare with the calculations done in R: plot(stand(tpa.eraint[,1]),lwd=3,col='grey') # From PMEL Ferret lines(stand(prtot.eraint[,1]),col='red',lty=3) # From R print(mean(tpa.eraint[,1])*0.5) print(mean(prtot.eraint[,1]))
The comparisons verified the consistency of the calulations in R and with Ferret, albeit with different units (not shown).
# Compare total precipitation between 50S-50N plot(stand(tpa.mon[,1]),lwd=3,col='grey') # From TRMM lines(stand(tpa.eraint[,1]),col='red',lty=3) # From ERAINT
Comparison for monthly precipitation area in ERAINT and TRMM (not shown)
# Compare precipitation area from monthly mean precipitation plot(stand(tpa.mon[,2]),lwd=3,col='grey') # From TRMM lines(stand(tpa.eraint[,2]),col='red',lty=3) # From ERAINT
CMIP 5 simulations
Estimate monthly aggregated precipitation area from the RCP4.5 CMIP5 experiment. These are only crude indices since the spatial resolution varies from model to model and the grid boxes represent an area average of the rainfall while the rain in reality only falls over a fraction of the area.
PareaCMIP() -> Parea.cmip5
Total precipitation
demo('global.t2m.gcm',ask=FALSE) ptot.cmip <- globalmean(path='~/CMIP5.monthly/rcp45', annual=TRUE,pattern='pr_',param='pr',anomaly=FALSE,FUN='sum') save(ptot.cmip,file='ptot.cmip.rda')
Loading the data for the analysis
For the analysis, the data/indices need to be loaded in memory:
data(Parea) data(Parea.month) # Only precipitation area data(Parea.eraint.month) # Only precipitation area data(prtot.eraint)
Some basic statistics such as the mean area of precipitation (km$^2$), wet-day mean precipitation (mm/day) and total precipitation (Ktons) estimated for the 50S-50N according to the TRMM data:
colMeans(Parea)
The analysis
Figures 1-4.
This code is provided for making the analysis more transparent and easier to replicate. However, the plotsa are shown in the main paper.
Fig 1 shows a map of a random day of precipirtation just to illustrate what is meant by the precipitation area (blue regions).
Fig1()
Fig 2 shows the daily precipitation area estimated from TRMM data.
Fig2()
t1 <- start(Parea) t2 <- end(Parea) A.p <- window(trend(Parea[,1]),start=t1,end=t2) ## Lower-Upper range print(range(Parea[,1])) ## P-value for trend analysis print(trend.pval(Parea[,1])) ## Trend: min and max print(range(A.p)) ## Change in percentage print(100*diff(range(A.p))/max(A.p)) ## Trend as percentage of the total area: print(round(100*range(A.p)/attr(Parea,'total area'),2))
Fig 3 compares monthly mean area estimated from daily area (TRMM and ERAINT) and monthly area estimated from monthly precipitation TRMM and ERAINT). Also compare monthly area between 50S--50N with global value. Fig 3 shows the daily precipitation area indices: both the daily indices and $A_P$ estimated from monthly mean precipitation for both. TRMM and ERAINT. The monthly aggregate differs from the daily area which seems to be shifted geographically from day to day, e.g. through storm tracks.
Fig3()
The monthly mean of daily indices do not follow the indeces of monthly mean precipitation. It is not expected that they should correspond, as they represent different aspects.
The time scale of preciptation is approximately hours, and often the precipitation does not move very far over the course of a day. Nevertheless, the 24-hr precipitation also involves moving storms and weather fronts to some extent, and are likely to over-estimate the area of precipitation.
The area of monthly mean precipitation is expected to be affected by the advection and migration of phenomena producing precipitation. For instance, the daily precipitation area $A_P$ could be constant, but the precipitation phenomena may move around so that the area of the monthly mean precipitation would appear to be larger. Even when setting a threshold $x_0$ for defining precipitation, the longer time scale will blur the area.
The spatial resolution too has an effect on the estimated area for a similar reason as the longer time scales. It is expected that it only rains in a fraction of the grid box, and the higher spatial resolution, the more precise we can expect the results to be.
The threshold $x_0$ was set to 1mm/day for calculation of the daily $A_P$but $x_0 = 4.5 mm/day$ for $A_P$ derived from monthly mean precipitation, except for ERANIT which had a coarser spatial resolution.
Fig 4 presents precipitation area indices for the CMIP5 RCP4.5 simulations were estimated from monthly mean precipitation $A_P = \sum_i a_i P_i\ \forall P_i > x_0$. The colour of the curves indicate their spatial resolution marked by the legend (number represents degree in logitude; red higher and blue lower spaital resolution).
Fig4()
Fig 4 also provides a summary of the trends indicated by the GCMs (%/decade). The mean increase over 2000--2100 is 2.4%.
Precipitation area indices for the CMIP5 RCP4.5 simulations were based on monthly mean precipitation. The CMIP5 results suggest an increase in the precipitation area for the monthly mean precipitation. Also the TRMM-based precipitation area based on monthly mean precipitation hinted towards an increase over time, suggesting more variable rainfall patterns given that the daily $A_P$ has declined. However, it is uncertain whether the GCMs also suggest a reduced $A_P$ for daily precipitation in the future.
For the CMIP5 data, $x_0 = 4 mm/day$ for calculating $A_P$.
Table GCM and mean precipitation area
Check how the mean precipitation area in (1000km)^2 varies wih the GCM and spatial resolution:
data(Parea.cmip5) y <- colMeans(Parea.cmip5)*1.e-6 names(y) <- attr(Parea.cmip5,'model_id') print(y)
Test to see if there is a systematic dependency on the model spatial resolution.
## ANOVA to analyse potential systematic relationship between model resolution and mean precipitation area data("IPCC.AR5.Table.9.A.1") nm1 <- tolower(gsub('.','',gsub('-','',IPCC.AR5.Table.9.A.1$Model.Name),fixed=TRUE)) nm2 <- tolower(gsub('.','',gsub('-','',attr(Parea.cmip5,'model_id')),fixed=TRUE)) ##x <- IPCC.AR5.Table.9.A.1$Horizontal.Grid[charmatch(nm2,nm1)] x <- cmipgcmresolution()[charmatch(nm2,nm1)] rbind(IPCC.AR5.Table.9.A.1$Model.Name,as.character(IPCC.AR5.Table.9.A.1$Horizontal.Grid),x) resbias <- data.frame(y=y,x=x) print(summary(lm(y ~ x, data=resbias)))
The results from the ANOVA suggests that there is a tendency for smaller mean precipitation area for models with higher resolution since the regression coefficient $x$ is negative. The relationship is statistically significant at the 5\%-level. One issue is that different spatial resolution may mean different values, as they reflect area mean values for the grid box. This may also influence the analysis as a constant threshold value $x_0$ for all these models may select precipitation events differently for models with high and low spatial resolution, unless the regridding of the models has taken the area-mean into account.
The resolution in this case was derived from a table of model details from the IPCC AR5, and the data here had all been regridded to the same spatial grid in the KNMI Climate Explorer https://climexp.knmi.nl/.
Supporting material and additional analysis
Rainfall area over ocean
The precipitation area over ocean was esimated by masking out all land/contentent grid boxes. This index is the third in the Parea object:
timeseries(Parea,is=2) print(trend.coef(Parea[,2])) print(trend.pval(Parea[,2])) print(range(window(trend(Parea[,2]),start=t1,end=t2)))
There was a clear decrease in $A_P$ over the oceans.
Rainfall area over land
Consequently, the rainfall area over land was estimated by masking all ocean grid points:
timeseries(Parea,is=3) print(trend.coef(Parea[,3])) print(trend.pval(Parea[,3])) print(range(window(trend(Parea[,3]),start=t1,end=t2)))
The total rain area over land is minor compared to the rain area over the oceans. Furthermore, the trend is much weaker and not statistically significant at the 5\%-level. There is a more pronounced annual cycle in $A_P$ over land than over the oceans, which may be related to monsoon-type systems.
Total precipitation from the TRMM data
The 50S-50N total precipitation amount was estimated by summing up all grid boxes with rainfall for each day:
timeseries(Parea,is=5,denominator=1e6,ylab='Gton') ## Lower-Upper range print(range(Parea[,5])/1.0e6) ## Check period print(c(start(Parea),end(Parea))) ## P-value for trend analysis print(trend.coef(Parea[,5])) print(trend.pval(Parea[,5])) ## Trend: min and max: ocean and land print(range(window(trend(Parea[,5]),start=t1,end=t2)/1e6)) ## Trend: min and max: ocean print(trend.coef(Parea[,6])) print(trend.pval(Parea[,6])) print(range(window(trend(Parea[,6]),start=t1,end=t2)/1e6)) ## Trend: min and max: land print(trend.coef(Parea[,7])) print(trend.pval(Parea[,7])) print(range(window(trend(Parea[,7]),start=t1,end=t2)/1e6))
The total daily volume of rainfall between 50S--50N exhibits a wide range of daily fluctuations between 868 -- 1494 Gton, but there is no strong long-term trend. An increasing trend is seen in the plot from 1120 -- 1150 Gton.
The total precipitation from ERAINT:
## ERAINT plot(prtot.eraint[,1]/1e9,ylab='Gton',main='ERAINT total precipitation 50S-50N',col='red',lwd=2) lines(trend(prtot.eraint[,1]/1e9)) ## Trend: min and max print(trend.coef(prtot.eraint[,1])) print(trend.pval(prtot.eraint[,1])) print(range(window(trend(prtot.eraint[,1]),start=t1,end=t2)/1e9))
Compare the total precipitation between 50S--50N with the rest of the globe:
print(colMeans(prtot.eraint)/1.0e9)
Rainfall amount over ocean
This index can be split into the total rainfall amount over the oceans by masking out the land grid boxes (not shown)
## The units are mm*km^2 -> 1000m^3 -> million kg or 1000 tons (1 Kton). timeseries(Parea,is=6,denominator=1e6,ylab='Gton')
The contribution to the total precipitation from the oceans exhibits a slightly more pronounced increasing trend than the total oceans+land amounts do.
The fraction of earth's area between 50S--50N is $\int$ $2 \pi a \cos (\phi) a d \phi/4 \pi a^2 = [\sin (\phi_2) - \sin (\phi_1)]/2$ where $\phi= \pm \pi*50/180$.
print(paste(round(100*sin(pi*50/180)),"% of earth's surface area is between 50S-50N"))
Rainfall amount over land
The corresponding index for the total land precipitation is (not shown)
timeseries(Parea,is=7,denominator=1.0e6,ylab='Gton')
The TRMM data suggests a declining trend in the total precipitation amount land in contrast to the precipitation over ocean and the sum oceans+land.
Precipitation intensity from TRMM
The rain intensity $\mu$ can be estimated by taking the average of all grid boxes with rain for each day:
timeseries(Parea,is=4,denominator=1) ## Trend: min and max print(trend.coef(Parea[,4])) print(trend.pval(Parea[,4])) print(round(range(window(trend(Parea[,4]),start=t1,end=t2)),2))
The graph reveals an increasing trend in $\overline{P}$, and the trend analysis is applied to a long series of daily values which makes it statistical significant at the 5\%-level (p-value $\approx 10^{-21}$).
A combination of a declining rainfall area over time with increasing precipitation amount is consistent with an increasing precipitation intensity: $P = \overline{P} * A_P$.
Check - covariation with the global mean temperature?
Assess the scaling relation between observed precipitation area and the global mean temperature through regression analysis ($A_P = \beta_0 + \beta_1 T + \eta$) where regression coefficient $\beta_1$ describes the slope (rate of change in $A_P$ with the global mean temperature) and $\beta_0$ the zero-intersect and $\eta$ being a nooise-term.
require(esd) hc4 <- HadCRUT4() ## The precipitation area: XY <- merge(annual(hc4),annual(Parea[,1]),all=FALSE) cal <- data.frame(x=coredata(XY[,1]),y=coredata(XY[,2])) pax <- lm(y ~x, data=cal) print(summary(pax)) ## Total precipitation: XY <- merge(annual(hc4),annual(Parea[,5]),all=FALSE) cal <- data.frame(x=coredata(XY[,1]),y=coredata(XY[,2])) print(summary(lm(y ~x, data=cal)))
An ANOVA on the TRMM annual mean values indicates strong covariance between the global mean temperature and the precipitation area, with $R^2$ of 0.6 and a negative regression coefficient $\beta_1 = -18 \times 10^6 km^2/$degree K. (statistically sign. at the 0.1\%-level). A similar analysis for the total precipitation suggests little covariance (not shown).
## Compare annual cycles xy <- merge(stand(window(hc4,start=t1,end=t2)),stand(Parea[,1]),all=TRUE) plot(aggregate(xy,month,mean,na.rm=TRUE),col=c('blue','red'),plot.type='single')
The there is little covariance between the mean seasonal cycles of HadCRUT4 and $P_A$.
Naive extrapolation for future daily $A_P$ based on global mean temperature change
data("global.t2m.cmip5") Pacmip5 <- global.t2m.cmip5 for (i in 1:dim(global.t2m.cmip5)[2]) { pre <- data.frame(x=coredata(global.t2m.cmip5[,i])) coredata(Pacmip5)[,i] <- predict(pax,newdata=pre) } plot(Pacmip5*1.0e-6,plot.type='single',main='Extrapolated precipitation area', ylab='million km^2',ylim=c(0,200),col=rgb(0,0,0,0.2)) grid() print(100*(mean(window(Pacmip5,start=2079,end=2099)) - mean(window(Pacmip5,start=1998,end=2016))) /mean(window(Pacmip5,start=1998,end=2016)))
A naive extrapolation assuming that (a) the regression analysis between the global mean temperature and precipitation area captures a true dependency and (b) that this dependency is stationary and holds in the future suggests a reduction in $A_P$ from near 110 million km$^2$ in the 20th century to 75 million km$^2$ by 2100. The CMIP5 ensemble mean suggests around a 28\% decrease in $A_P$.
Total atmospheric water content from ERAINT
The total 50S-50N atmospheric water content was estimated from the ERAINT reanalysis
data(H2O) plot(H2O[,1]/1e9,ylab='Gton') ## The 50S-50N water content is much higher than for the high latitude regions ## "kg km**-2" agregated (sum) over the area lines(trend(H2O[,1]/1e9)) grid() ## Trend: min and max print(trend.coef(H2O[,1])) print(trend.pval(H2O[,1])) print(range(window(trend(H2O[,1]),start=t1,end=t2)/1e9)) print(c(start(H2O),end(H2O))) colMeans(H2O)/1e9 ## Gtons
The total atmospheric water exhibits weak long-term changes according to the ERAINT reanalysis and shows a hint of an increase over time.
The surface moisture flux $Q_s$ from ERAINT is the same as evaporation $E_s$ and the declining trends (negative numbers) indicates an increased evaporation over time. This is consistent with increased $P_{tot}$ in the TRMM data. Negative values indicate flux out of the oceans.
data(Qs) plot(annual(Qs,FUN='mean'),ylab='Gton') ## The 50S-50N moisture flux is much higher than for the high latitude regions ## Units are mm/day * km^2 which is equivalent to 1000m^3/day or kilo-tons/day lines(trend(annual(Qs,FUN='mean'))) grid() ## Giga-tons/day - divide by 1e6 ## Trend: min and max print(trend.coef(Qs[,1])) print(trend.pval(Qs[,1])) print(range(window(trend(Qs[,1]),start=t1,end=t2))) print(c(start(Qs),end(Qs))) mean(Qs)
Table on the amount of evaporation and precipitation - water fluxes
Table 1:
print('Moisture flux (ERAINT; Gton/day)') print(round(mean(Qs))) print('Column water (ERAINT; Gton)') print(round(colMeans(H2O)*1.e-9)) print('Precipitation (ERAINT; Gton)') print(round(colMeans(prtot.eraint)*1.e-9))
print('Mean total precip area 50S-50N (Gton/day)') print(round(mean(Parea[,1]*1.e-6))) print('Mean total precip area over ocean 50S-50N (Gton/day)') print(round(mean(Parea[,2]*1.e-6))) print('Mean total precip area 50S-50N over land (Gton/day)') print(round(mean(Parea[,3]*1.e-6))) print('Mean total precip amount 50S-50N (Gton/day)') print(round(mean(Parea[,5]*1.e-6))) print('Mean total precip amount over ocean 50S-50N (Gton/day)') print(round(mean(Parea[,6]*1.e-6))) print('Mean total precip amount 50S-50N over land (Gton/day)') print(round(mean(Parea[,7]*1.e-6)))
Table 2:
h2o <- window(trend(H2O[,1]),start=t1,end=t2)/1e15 print(round(c(h2o[1],h2o[length(h2o)]))) qs <- window(trend(Qs),start=t1,end=t2) print(round(c(qs[1],qs[length(qs)]))) y <- window(trend(Parea[,1]),start=t1,end=t2)/1e6 print(round(c(y[1],y[length(y)]))) y <- 100*window(trend(Parea[,1]),start=t1,end=t2)/attr(Parea,'total area') print(round(c(y[1],y[length(y)]))) y <- window(trend(Parea[,2]),start=t1,end=t2)/1e6 print(round(c(y[1],y[length(y)]))) y <- window(trend(Parea[,3]),start=t1,end=t2)/1e6 print(round(c(y[1],y[length(y)]))) y <- window(trend(Parea[,5]),start=t1,end=t2)/1e6 print(round(c(y[1],y[length(y)]))) y <- window(trend(prtot.eraint[,1]),start=t1,end=t2)/1e9 print(round(c(y[1],y[length(y)]))) y <- window(trend(Parea[,6]),start=t1,end=t2)/1e6 print(round(c(y[1],y[length(y)]))) y <- window(trend(Parea[,7]),start=t1,end=t2)/1e6 print(round(c(y[1],y[length(y)]))) y <- window(trend(Parea[,4]),start=t1,end=t2) print(round(c(y[1],y[length(y)]),1)) y <- window(trend(Parea.month[,1])/1e6,start=t1,end=t2) print(round(c(y[1],y[length(y)]),1)) y <- window(trend(tpa.eraint[,3]),start=t1,end=t2)/1e12 # Average of two print(round(c(y[1],y[length(y)]),1)) index(Parea.eraint.month) <- as.Date(index(Parea.eraint.month)) y <- window(trend(Parea.eraint.month),start=t1,end=t2)/1e6 print(round(c(y[1],y[length(y)]),1))
Previous estimates of global precipitation flux
Table 1 in Zektser et al (1993) http://www.sciencedirect.com/science/article/pii/0022169493901829
## Example in the text: evaporation print(294*128) ## km^3/year - should be 38,000 km^3 according to the article ## As kg/day: 1 km^3 = 1.oe9 tons = 1 Gton. print(294*128/365.25) ## Precipitation over land print(paste('Precipitation: 834 mm/year=',round(834*128/365.25),'Gton/day'))
Estimate of the global mean precipitation from Kiel and Trenberth: 2.46 -- 2.90 mm/day http://journals.ametsoc.org/doi/pdf/10.1175/1520-0477%281997%29078%3C0197%3AEAGMEB%3E2.0.CO%3B2
print(paste('2.46mm/day =',round(2.46*1.0e-3 * 4*pi*6378000^2*1.0e-9),'Gton/day')) print(paste('2.46mm/day =',round(2.90*1.0e-3 * 4*pi*6378000^2*1.0e-9),'Gton/day'))
Legates (1995): 2.6 to 3.1 mm/day http://onlinelibrary.wiley.com/doi/10.1002/joc.3370150302/full
print(paste('2.6mm/day =',round(2.6*1.0e-3 * 4*pi*6378000^2*1.0e-9),'Gton/day')) print(paste('3.1mm/day =',round(3.1*1.0e-3 * 4*pi*6378000^2*1.0e-9),'Gton/day'))
Test the correlation between the precipitation area and galactic cosmic rays
Compare the cosmic rays from Climax with the precipitation area (not shown)
require(replicationDemos) data(gcr) t <- paste(substr(gcr$Date,7,10),substr(gcr$Date,1,2),substr(gcr$Date,4,5),sep='-') climax <- zoo(as.numeric(as.character(gcr$Climax)),order.by=as.Date(t)) plot(climax,Parea[,1]) gcrpa <- merge(climax,Parea[,1],all=FALSE) ok <- is.finite(gcrpa[,1]) & is.finite(gcrpa[,2]) print(cor.test(gcrpa[ok,1],gcrpa[ok,2]))
Sanity check - test of consistency
Compare the total precipitation from the daily data and the monthly data
P.tot <- merge(aggregate(Parea[,5],year,mean), aggregate(prtot.eraint[,1]/1000,year,mean),all=TRUE) plot(P.tot/1.0e6,plot.type='single',col=c('black','red'),led=c(2,2))
Estimate of the portion of the vertical energy flow is due to convection from the total precipitation
The precipitation was assumed to be a by-product of moist convection with an energy content equivalent to the release of latent heat of evaporation $P L_E$:
$f_{conv}= \frac{P L_e}{\pi a^2 (1-A) S_0}$
a <- 6378000 # Earth's area (m) A = 0.3 # The albedo S0 <- 1361 # The Solar constant (W/m²) Le <- 2264.76 * 1000 # Latent heat of evaporation (https://en.wikipedia.org/wiki/Latent_heat) J/kg tot.E <- pi * a^2 * (1-A) * S0 * 24 * 3600 # Energy over 24 hrs P <- sum(colMeans(prtot.eraint))*1000 frac.E <- P*Le/tot.E * 100 print(paste('Convection is responsible for ',round(frac.E,2),'percent of the vertical heat flow'))
Trend anlysis of the proportn of vertical energy flow being connected to convection
fE <- (prtot.eraint[,1]+prtot.eraint[,2]+prtot.eraint[,3])*Le*100000/tot.E plot(fE) grid() lines(trend(fE)) prop <- window(trend(fE),start=t1,end=t2) print(round(c(prop[1],prop[length(prop)]),1))
For the TRMM data:
y <- window(trend(Parea[,5]),start=t1,end=t2) print(paste('Percentage energy flow associated with convection',round(c(y[1],y[length(y)])*Le*1e6/tot.E*100),'%'))
Estimate the globally summed precipitation from CMIP5:
data(ptot.cmip) plot(zoo(ptot.cmip)/1e6,plot.type='single',ylab='Gton') pclim <- colMeans(window(ptot.cmip,start=1980,end=2009))/1e6 # mean total precipitation print(summary(pclim*Le/tot.E * 100)) print(summary(pclim)) # In gigatons Ptrends <- apply(window(ptot.cmip,start=2000,end=2100)/1e6,2,trend.coef) print(summary(Ptrends)) Ftrends <- apply(100*window(ptot.cmip,start=2000,end=2100)*Le*1e6/tot.E,2,trend.coef) ## Change in percent of energy flow connected with convection between 2000 and 2100 print(summary(Ftrends*10))
Acknowledgements:
The PMEL/NOAA software ferret http://ferret.pmel.noaa.gov/Ferret/ and David Pierce's ncview were used to check and verify the netCDF files and the calculations. The function globalsumand ferret indicated similar results for the monthly mean precipitation totals.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.