Linde Temperature / Water balance analysis

Use daily climate station time series from DWD, downloaded here: https://opendata.dwd.de/climate_environment/CDC/observations_germany/climate/daily/kl/historical/

Chose stations - No 6265 Wusteritz (01.5.2004 - current) - Data missing - start in 2005, same period as Muehlhausen is available

filename = "tageswerte_KL_06265_20040501_20221231_hist/produkt_klima_tag_20040501_20221231_06265.txt" # Brandenburg/Havel

There were 23 warnings (use warnings() to see them)

Check which period is covered, whether years are complete

data$date <- as.Date(as.character(data$MESS_DATUM),format='%Y%m%d')
years <- as.numeric(unique(format(data$date,'%Y')))

data$year <- format(data$date,'%Y')
data$month <- format(data$date,'%m')
data$daymonth <- format(data$date,'%d%m')
data$day <- format(data$date,'%d')

(days_per_year <- rle(data$year))

Run Length Encoding
  lengths: int [1:19] 245 365 365 365 366 365 365 365 366 365 ...
  values : chr [1:19] "2004" "2005" "2006" "2007" "2008" "2009" "2010" "2011" "2012" "2013" "2014" "2015" "2016" ...

Select only a subset of years: 2005 - 2021 (will yields 16 hydrological years) which is aligned with the data availability of Hainich. Incidentally, it covers 2014-2021 (8 years), and the preceding 8 years.

data = subset(data, ((year >= 2005) & (year<= 2021)))
years <- as.numeric(unique(format(data$date,'%Y')))
(nyears = length(years))

[1] 17

days_per_year <- rle(data$year)
days_per_year[[1]]

 [1] 365 365 365 366 365 365 365 366 365 365 365 366 365 365 365 366 365

remove Feb 29

ind_feb29<-which(as.character(data$date,'%d%m') == c("2902"))
data<-data[-ind_feb29,]

Check whether there are any missing data that need to be taken care of.

(a <- which(data$TMK == -999)) # daily mean temp, 

integer(0)

(a <- which(data$TXK == -999)) # daily max temp, 

integer(0)

(a <- which(data$TNK == -999)) # daily min temp, 

integer(0)

(data$date[a]) # need to check surrounging stations

Date of length 0

(b <- which(data$RSK == -999)) # daily precipitation

[1] 3876 3877 3878 3879 3880

(data$date[b]) # need to check surrounging stations

[1] "2015-08-14" "2015-08-15" "2015-08-16" "2015-08-17" "2015-08-18"

data$RSK[b] = c(4.7,0,24.5,0.3,0.5)
(b <- which(data$RSK == -999)) # daily precipitation

integer(0)

Data for precipition are missing 5 consecutive days in 2015. Replace with values from teh station near the Linde site: Berge

STATIONS_ID MESS_DATUM QN_3 FX FM QN_4 RSK RSKF SDK SHK_TAG NM VPM PM TMK UPM TXK TNK TGK eor date year month daymonth day

   5825   20150810   10  9.1 2.2   10  0.0    0 11.6    -999 -999 17.0 -999 24.2  57 32.6 14.3 13.8 eor 2015-08-10 2015    08     1008  10

   5825   20150811   10 11.1 3.4   10  0.0    0 12.7    -999 -999 19.4 -999 24.2  66 31.7 18.2 17.4 eor 2015-08-11 2015    08     1108  11

   5825   20150812   10  9.2 3.4   10  0.1    4  5.3    -999 -999 18.3 -999 21.2  74 27.8 17.1 16.5 eor 2015-08-12 2015    08     1208  12

   5825   20150813   10  8.1 2.9   10  0.0    0 10.5    -999 -999 15.9 -999 22.7  60 28.7 17.4 16.4 eor 2015-08-13 2015    08     1308  13

** 5825 20150814 10 11.2 3.8 10 4.7 4 12.0 -999 -999 16.6 -999 26.2 51 34.0 19.8 18.6 eor 2015-08-14 2015 08 1408 14

   5825   20150815   10 12.5 2.4   10  0.0    0  7.1    -999 -999 21.6 -999 22.3  81 29.1 18.6 17.5 eor 2015-08-15 2015    08     1508  15

   5825   20150816   10 10.9 2.3   10 24.5    4  6.3    -999 -999 20.5 -999 21.8  82 30.0 17.6 17.1 eor 2015-08-16 2015    08     1608  16

   5825   20150817   10 10.4 2.5   10  0.3    4  1.2    -999 -999 18.2 -999 22.1  70 26.8 19.0 17.9 eor 2015-08-17 2015    08     1708  17

   5825   20150818   10  5.6 2.2   10  0.5    4  0.0    -999 -999 18.1 -999 17.7  89 20.2 16.3 16.1 eor 2015-08-18 2015    08     1808  18 **

   5825   20150819   10  7.3 2.5   10  0.0    0  5.0    -999 -999 17.1 -999 16.6  91 20.2 13.9 12.5 eor 2015-08-19 2015    08     1908  19

   5825   20150820   10  6.7 1.5   10  0.0    0 11.2    -999 -999 12.8 -999 18.1  67 25.0 12.0 10.2 eor 2015-08-20 2015    08     2008  20

Add a column indicating the daily difference between minimum and maximum temperature. This is a measure for the potential evapotranspiration.

data$DeltaT <- data$TXK - data$TMK

Calculation of the hydrological variables

Potential ET

Calculation, taking the equation of Hargreaves-Samani (1985), see review by Hargreaves and Allen (2003). This takes a bit.

library(sirad)

a <- extrat(dayOfYear(data$date), Lat_data) # calculate extraterrestrial radiation for the specific day and latitude
data$Ra <- a[[1]]
data$lEpot <- 0.0023 * data$Ra * (data$TMK + 17.8) * sqrt(data$DeltaT) # calculate potential latent heat flux in MJ/m2d
data$Epot <- (data$lEpot/2.5e6)*1e6# calculate potential ET in mm/d assume latent heat of evaporation l = 2.5e6 J/kg and density of water of 1000 kg/m3

Atmospheric water deficit

The water deficit can either be calculated over each the hydrological year (October - September) separately or over the entire period.

over each hydrological year separately

hydro.startyear = years[1]
hydro.endyear = years[length(years)]

hydro.startdate = as.Date(paste(hydro.startyear,"10","01", sep="-"), "%Y-%m-%d")
hydro.enddate = as.Date(paste(hydro.endyear,"09","30", sep="-"), "%Y-%m-%d")

#hydro.startdate = as.POSIXlt(paste(hydro.startyear,"10","01", sep="-"), "%Y-%m-%d")
#hydro.enddate = as.POSIXlt(paste(hydro.endyear,"09","30", sep="-"), "%Y-%m-%d")

h.data <- subset.data.frame(data, (date >= hydro.startdate & date <= hydro.enddate))

h.data$def <- NA


for (year in years[1:length(years)-1]) {
        startday = as.Date(paste(year,"10","01", sep="-"), "%Y-%m-%d")
        endday = as.Date(paste(year+1,"09","30", sep="-"), "%Y-%m-%d")
        ind.thisyear <- which(h.data$date >= startday & h.data$date <= endday)
        h.data$def[ind.thisyear] <- cumsum(h.data$RSK[ind.thisyear]-h.data$Epot[ind.thisyear]) 

}

Basic features of the temperature dataset

Check whether the distribution of rainfall changed between

• 2014-2017 and

• 2018-2021.

Since the observations from Linde cover a shorter period (2020-2021), create an alterative subset, to be able to look only at some of the analysis for this smaller range.

• 2020 and

For orientation plot the dayliy climatology over the entire timeseries, 2005-2021

# subset for periods
data_2014_2017 <- data[which(data$year>=2014 & data$year<=2017),] 
data_2018_2021 <- data[which(data$year>=2018 & data$year<=2021),] 

data_2020 <- data[which(data$year==2020),] 
data_2021 <- data[which(data$year==2021),] 

# climatology over the entire time series
T.clim <- aggregate(TMK ~ day+month, data = data, FUN = mean)
plot(T.clim$TMK, ylim = c(-7, 28),
     xlab = "Doy of year",
     ylab = "Average daily temperature in deg C")

# averge over periods
T.clim.2014_2017 <- aggregate(TMK ~ day+month, data = data_2014_2017, FUN = mean)
points(T.clim.2014_2017$TMK, col="blue")

# 2018-2021
T.clim.2018_2021 <-aggregate(TMK ~ day+month, data = data_2018_2021, FUN = mean)
points(T.clim.2018_2021$TMK, col="red")
legend(x = "topleft",          # Position
       legend = c("clim", "2014-2017", "2018-2021"),  # Legend texts
       pch = c(1,1,1), # point type
       col = c("black", "blue", "red")) # point colors

plot(T.clim$TMK-T.clim.2014_2017$TMK, col="blue",
     xlab=c("Day of year"),
     ylab = "Temperature anomaly in degC")
points(T.clim$TMK-T.clim.2018_2021$TMK, col = "red")

legend(x = "topleft",          # Position
       legend = c("2014-2017", "2018-2021"),  # Legend texts
       pch = c(1,1), # point type
       col = c("blue", "red")) # point colors



par(mfrow=c(1,2))

bins = seq(-6,6,1)
hist(T.clim.2014_2017$TMK-T.clim$TMK,
     xlab = "Deviation in temperature from climatology",
     main = "2014-2018",
     breaks = bins)
hist(T.clim.2018_2021$TMK-T.clim$TMK,
     xlab = "Deviation in temperature from climatology",
     main = "2018-2021",
     breaks = bins)

(quantile(T.clim.2014_2017$TMK-T.clim$TMK, c(0.1,0.5,0.9)))

       10%        50%        90% 
-1.5402941  0.2264706  2.1329412 

(quantile(T.clim.2018_2021$TMK-T.clim$TMK, c(0.1,0.5,0.9)))

       10%        50%        90% 
-0.9967647  0.6132353  2.8008824 

# temperature anomalies
data$Tclim <- NA
ind <- which(data$year>years[1])
data$Tclim[ind] <- rep(T.clim$TMK,length(ind)/365) 
data$Tano <- data$TMK - data$Tclim

plot(data$date,data$Tano,
     xlab = "Date",
     ylab = "Temperature anomalies in degC")

boxplot(Tano ~ year, data=data,
        xlab = "Year",
        ylab = "Temperature anomalies in degC")

lm.Tanom <- lm(Tano ~ as.numeric(year), data = data)
summary(lm.Tanom)

Call:
lm(formula = Tano ~ as.numeric(year), data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.9483  -2.4361  -0.0775   2.4687  10.2643 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -132.79228   20.84488  -6.370 2.03e-10 ***
as.numeric(year)    0.06597    0.01035   6.372 2.01e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.647 on 5838 degrees of freedom
  (365 observations deleted due to missingness)
Multiple R-squared:  0.006907,  Adjusted R-squared:  0.006737 
F-statistic: 40.61 on 1 and 5838 DF,  p-value: 2.006e-10

Both periods were warmer than the entire reference periods from 2005-2021, while 2018-2021 was overall roughly 0.05 degree warmer than the previous 2014-2017 period. There are no strong effects in temperature per se. ALso the anomalies are slightly increasing with time by 0.066 degree per year.

Difference in the daily temperature indicating cloudiness.

(delta_T = T_max - T_min)

DeltaT.clim <- aggregate(DeltaT ~ day+month, data = data, FUN = mean)
plot(DeltaT.clim$DeltaT, ylim = c(1,10))

# averge over periods
DeltaT.clim.2014_2017 <- aggregate(DeltaT ~ day+month, data = data_2014_2017, FUN = mean)
points(DeltaT.clim.2014_2017$DeltaT, col="blue")

# averge over periods
DeltaT.clim.2018_2021 <- aggregate(DeltaT ~ day+month, data = data_2018_2021, FUN = mean)
points(DeltaT.clim.2018_2021$DeltaT, col="red")

DeltaT.clim.2020 <- aggregate(DeltaT ~ day+month, data = data_2020, FUN = mean)
DeltaT.clim.2021 <- aggregate(DeltaT ~ day+month, data = data_2021, FUN = mean)


bins = seq(0.1,12,0.25)

h.DeltaT_daily.clim <- hist(DeltaT.clim$DeltaT, breaks = bins, plot = FALSE)
h.DeltaT_daily.2014_2017 <- hist(DeltaT.clim.2014_2017$DeltaT, breaks = bins, plot = FALSE)
h.DeltaT_daily.2018_2021 <- hist(DeltaT.clim.2018_2021$DeltaT, breaks = bins, plot = FALSE)
h.DeltaT_daily.2020 <- hist(DeltaT.clim.2020$DeltaT, breaks = bins, plot = FALSE)
h.DeltaT_daily.2021 <- hist(DeltaT.clim.2021$DeltaT, breaks = bins, plot = FALSE)

plot(h.DeltaT_daily.clim$mids, cumsum(h.DeltaT_daily.clim$density)/4, 
     xlab = "Daily temperature difference, degC",
     ylab = "cumulative frequency")
points(h.DeltaT_daily.2014_2017$mids, cumsum(h.DeltaT_daily.2014_2017$density)/4, 
       pch = 19, col = "blue")

points(h.DeltaT_daily.2018_2021$mids, cumsum(h.DeltaT_daily.2018_2021$density)/4, 
       pch = 19, col = "red")
legend(x = "bottomright",          # Position
       legend = c("2005-2021", "2014-2017", "2018-2021"),  # Legend texts
       pch = c(1, 19, 19), # point types
       col = c("black", "blue", "red"), # point colors
       )           

plot(h.DeltaT_daily.clim$mids, cumsum(h.DeltaT_daily.clim$density)/4, 
     xlab = "Daily temperature difference, degC",
     ylab = "cumulative frequency")
points(h.DeltaT_daily.2020$mids, cumsum(h.DeltaT_daily.2020$density)/4, 
       pch = 19, col = "blue")

points(h.DeltaT_daily.2021$mids, cumsum(h.DeltaT_daily.2021$density)/4, 
       pch = 19, col = "red")
legend(x = "bottomright",          # Position
       legend = c("2005-2021", "2020", "2021"),  # Legend texts
       pch = c(1, 19, 19), # point types
       col = c("black", "blue", "red"), # point colors
       )           

deltaT <- as.data.frame(cbind(h.DeltaT_daily.clim$mids, cumsum(h.DeltaT_daily.clim$density)/4, cumsum(h.DeltaT_daily.2014_2017$density)/4, cumsum(h.DeltaT_daily.2018_2021$density)/4))
colnames(deltaT) <- c("delta_T", "freq_2005_2021", "freq_2014_2017", "freq_2018_2021")

write.csv(deltaT, file="daltaT_15yrs_6265.csv")

#hist(DeltaT.clim.2014_2017$DeltaT-DeltaT.clim$DeltaT)

#hist(DeltaT.clim.2018_2021$DeltaT-DeltaT.clim$DeltaT)

The difference between minumum and maxumim temperature increases in both periods, esply in times were temperature differences are high (in summer). Thus in general the variation in temperature at a daily scale has increased, indicating an increase in sunshine hours and decrease in cloudiness.

Potential Evapotranspiration, daily

Epot.clim <- aggregate(Epot ~ day+month, data = data, FUN = mean)
plot(Epot.clim$Epot, ylim = c(0,5),
     xlab = "Day of year",
     ylab = "Potential evapotranspiration in mm/day")
legend(x = "topright",          # Position
       legend = c("clim", "2014-2017", "2018-2021"),  # Legend texts
        pch = c(1, 1, 1), # point types
        col = c("black", "blue", "red")) # point colors

# averge over periods
Epot.clim.2014_2017 <- aggregate(Epot ~ day+month, data = data_2014_2017, FUN = mean)
points(Epot.clim.2014_2017$Epot, col="blue")

# averge over periods
Epot.clim.2018_2021 <- aggregate(Epot ~ day+month, data = data_2018_2021, FUN = mean)
points(Epot.clim.2018_2021$Epot, col="red")

bins = seq(0,5,0.1)

h.Epot_daily.clim <- hist(Epot.clim$Epot, breaks = bins, plot = FALSE)
h.Epot_daily.2014_2017 <- hist(Epot.clim.2014_2017$Epot, breaks = bins, plot = FALSE)
h.Epot_daily.2018_2021 <- hist(Epot.clim.2018_2021$Epot, breaks = bins, plot = FALSE)

plot(h.Epot_daily.clim$mids, cumsum(h.Epot_daily.clim$density)/10, 
     xlab = "Daily potential evapotranspiration, mm/d",
     ylab = "cumulative frequency")
points(h.Epot_daily.2014_2017$mids, cumsum(h.Epot_daily.2014_2017$density)/10, 
       pch = 19, col = "blue")

points(h.Epot_daily.2018_2021$mids, cumsum(h.Epot_daily.2018_2021$density)/10, 
       pch = 19, col = "red")
legend(x = "bottomright",          # Position
       legend = c("clim", "2014-2017", "2018-2021"),  # Legend texts
       pch = c(1, 19, 19), # point types
       col = c("black", "blue", "red"), # point colors
       )           

At the daily scale, the potential evapotranspiration has slightly increased especially in summer.

Potential evapotranspiration, precipitation weekly

library(tidyr)
# Aggregate to weekly data as 5 daily averages perform better than daily 
data$date <- as.POSIXlt(data$date)
data.w.Epot <- aggregate(data$Epot ~ format(date,"%Y-%W"), data=data, FUN = mean) # weekly data
names(data.w.Epot)[1:2] = c("date", "Epot")
data.w.Epot<-data.w.Epot %>% separate(date, c("year","week")) # tidyr

# subset for periods once more
data.w.Epot_2014_2017 <- data.w.Epot[which(data.w.Epot$year>=2014 & data.w.Epot$year<=2017),] 
data.w.Epot_2018_2021 <- data.w.Epot[which(data.w.Epot$year>=2018 & data.w.Epot$year<=2021),] 

Epot.clim <- aggregate(Epot ~ week, data = data.w.Epot, FUN = mean) # weekly climatology

Epot.clim.2014_2017 <- aggregate(Epot ~ week, data = data.w.Epot_2014_2017, FUN = mean) # weekly climatology 2014-2017
Epot.clim.2018_2021 <- aggregate(Epot ~ week, data = data.w.Epot_2018_2021, FUN = mean) # weekly climatology 2018-2021

plot(Epot.clim$Epot,
     xlab = "Calendar week",
     ylab = "Average PET in mm/d")
points(Epot.clim.2014_2017$Epot, 
       pch = 19, col = "blue")

points(Epot.clim.2018_2021$Epot, 
       pch = 19, col = "red")
legend(x = "topleft",          # Position
       legend = c("clim", "2014-2017", "2018-2021"),  # Legend texts
       pch = c(1, 19, 19), # point types
       col = c("black", "blue", "red"), # point colors
       )           

plot(cumsum(Epot.clim$Epot*7),
     xlab = "Calendar week",
     ylab = "Cumulated PET in mm")
points(cumsum(Epot.clim.2014_2017$Epot*7), 
       pch = 19, col = "blue")

points(cumsum(Epot.clim.2018_2021$Epot*7), 
       pch = 19, col = "red")
legend(x = "bottomright",          # Position
       legend = c("clim", "2014-2017", "2018-2021"),  # Legend texts
       pch = c(1, 19, 19), # point types
       col = c("black", "blue", "red"), # point colors
       ) 

NA
NA
NA

This gives the cumulative potential evapotranspiration, which has slightly increased in the second period compared to the first. The first is similar to the climatology.

Monthly precipitation and potential evapotranspiration

prec.monthly <- aggregate(RSK ~ month+year, data = data, FUN = sum)
clim.prec.monthly <- aggregate(RSK ~ month, data = prec.monthly, FUN = mean)
plot(clim.prec.monthly,
     ylim = c(0,150),
     type = "l",
     ylab = "Monthly flux, mm")

# period 2014-2017
prec.monthly.2014_2017 <- aggregate(RSK ~ month+year, data = data_2014_2017, FUN = sum)
clim.prec.monthly.2014_2017 <- aggregate(RSK ~ month, data = prec.monthly.2014_2017, FUN = mean)
lines(clim.prec.monthly.2014_2017,
     type = "l",
     col = "blue")

# period 2018-2021
prec.monthly.2018_2021 <- aggregate(RSK ~ month+year, data = data_2018_2021, FUN = sum)
clim.prec.monthly.2018_2021 <- aggregate(RSK ~ month, data = prec.monthly.2018_2021, FUN = mean)
lines(clim.prec.monthly.2018_2021,
     type = "l",
     col = "red")


# all years
Epot.monthly <- aggregate(Epot ~ month+year, data = data, FUN = sum)
clim.Epot.monthly <- aggregate(Epot ~ month, data = Epot.monthly, FUN = mean)

# period 2014-2017
Epot.monthly.2014_2017 <- aggregate(Epot ~ month+year, data = data_2014_2017, FUN = sum)
clim.Epot.monthly.2014_2017 <- aggregate(Epot ~ month, data = Epot.monthly.2014_2017, FUN = mean)

# period 2018-2021
Epot.monthly.2018_2021 <- aggregate(Epot ~ month+year, data = data_2018_2021, FUN = sum)
clim.Epot.monthly.2018_2021 <- aggregate(Epot ~ month, data = Epot.monthly.2018_2021, FUN = mean)

lines(clim.Epot.monthly,
     type = "l",
     lty = "dashed")

lines(clim.Epot.monthly.2014_2017,
     type = "l",
     col="blue",
     lty = "dashed")
lines(clim.Epot.monthly.2018_2021,
     type = "l",
     col="red",
     lty = "dashed")

legend(x = "topleft",          # Position
       legend = c("precipitation","Epot", "all", "2014-2017", "2018-2021"),  # Legend texts
       lty = c("solid", "dashed", "solid", "solid","solid"), # point types
       col = c("black", "black", "black", "blue", "red"), # point colors
       ) 

Precipitation varies more between the focus periods than potential evaporation.

Water balance, per hydrological year

water deficit = precipitation - potential evaporation

plot(h.data$date, h.data$def, 
#     ylim = c(-200,250),
     xlab = "Date",
     ylab = "P - PET since Oct 1 each year")

In this plot, all hydrological years are considered separately. The water deficit is set to zero at the beginning of each hydrological year (corresponds to the end of the vegetation period) at Oct 1.

We see that the variation within the years increased with time.

def.yearly.mat <- matrix(h.data$def, nrow=365, ncol=length(years)-1 )
h.years = years[-1]

def.min.year <- apply(def.yearly.mat,2,min)
def.max.year <- apply(def.yearly.mat,2,max)
def.range.year <- def.max.year - def.min.year
def.sd.year <- apply(def.yearly.mat,2,sd)  
def.mean.year <- apply(def.yearly.mat,2,mean)


plot(years[1:length(years)-1], def.range.year,
     ylab=c("Range of water availability"),
     xlab=c("Calendar year"))    

ind.2014_2021 = c(10:17)

lm.def.range <- lm(def.range.year ~ years[2:length(years)])
(summary(lm.def.range))

Call:
lm(formula = def.range.year ~ years[2:length(years)])

Residuals:
   Min     1Q Median     3Q    Max 
-98.94 -31.92 -10.66  17.77 141.80 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)
(Intercept)            -1375.2515  7037.1910  -0.195    0.848
years[2:length(years)]     0.8024     3.4950   0.230    0.822

Residual standard error: 64.44 on 14 degrees of freedom
Multiple R-squared:  0.003751,  Adjusted R-squared:  -0.06741 
F-statistic: 0.05271 on 1 and 14 DF,  p-value: 0.8217

lm.def.sd <- lm(def.sd.year ~ years[2:length(years)])
(summary(lm.def.sd))

Call:
lm(formula = def.sd.year ~ years[2:length(years)])

Residuals:
    Min      1Q  Median      3Q     Max 
-31.721  -6.586  -4.509   4.833  49.811 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)
(Intercept)            -783.5880  2302.9587  -0.340    0.739
years[2:length(years)]    0.4226     1.1438   0.369    0.717

Residual standard error: 21.09 on 14 degrees of freedom
Multiple R-squared:  0.009657,  Adjusted R-squared:  -0.06108 
F-statistic: 0.1365 on 1 and 14 DF,  p-value: 0.7173

#plot(years[1:length(years)-1], def.sd.year,
#     ylab="variation of water availability",
#     xlab="calendar year")    

#plot(years[1:length(years)-1], def.mean.year,
#     ylab=c("mean annual water deficit"),
#     xlab=c("calendar year"))    

#plot(years[1:length(years)-1], def.min.year,
#     ylab=c("minimum water availability"),
#     xlab=c("calendar year"))


#plot(years[1:length(years)-1], def.max.year, 
#     ylab=c("maximum water availability"),
#     xlab=c("calendar year"))

This plot shows the difference between the yearly maximum surplus and minimum water deficit, called “range in water deficit” (positive values of water deficit are surplus). The range decreases significantly in the recent years.

clim.def <- apply(def.yearly.mat,1,mean)

c2014_2017 <- which(h.years>=2014 & h.years<=2017)
c2018_2021 <- which(h.years>=2018 & h.years<=2021)

h.clim.date = h.data$date[1:365]
clim.def.mean.2014_2017 <- apply(def.yearly.mat[,c2014_2017],1,mean)
clim.def.min.2014_2017 <- apply(def.yearly.mat[,c2014_2017],1,min)
clim.def.max.2014_2017 <- apply(def.yearly.mat[,c2014_2017],1,max)

clim.def.mean.2018_2021 <- apply(def.yearly.mat[,c2018_2021],1,mean)
clim.def.max.2018_2021 <- apply(def.yearly.mat[,c2018_2021],1,max)
clim.def.min.2018_2021 <- apply(def.yearly.mat[,c2018_2021],1,min)


plot(h.clim.date, clim.def,
     ylim = c(-150, 200),
     xlab = "Date",
     ylab = "Water deficit/surplus = P-PET in mm")
points(h.clim.date, clim.def.mean.2014_2017,
       col="blue")

points(h.clim.date, clim.def.mean.2018_2021,
       col="red")
legend(x = "topleft",          # Position
       legend = c("clim", "2014-2017", "2018-2021"),  # Legend texts
       pch = c(1, 19, 19), # point types
       col = c("black", "blue", "red"), # point colors
       )

NA
NA

When aggreagating over the two periods, we see the general dynamics. There is a surplus of water in the winter, because ET is low, and a slight decrease in summer, when ET is larger, but also P increases. The averge multiannual surplus is 264 mm/year! The most striking difference between the first and second period is that the met water balance in 2014-2017 was above average (353 mm).

#plot(h.clim.date, clim.def,
#     ylim = c(-100, 150),
#     xlab = "date",
#     ylab = "water deficit in P-PET in mm",
#     type = "l")
#polygon(c(h.clim.date), rev(h.clim.date), 
#          c(clim.def.max.2014_2017, rev(clim.def.min.2014_2017)),
#        col = "#6BD7AF")

days = c(1:365)
daysmonth = c(31,30,31,31,28,31,30,31,30,31,31,30)
months = c("Oct","Nov","Dec","Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep")
axisticks = c(31/2,cumsum(daysmonth[1:11]) + daysmonth[2:12]/2)

plot(days, clim.def,
    xaxt = "n" ,
    ylim = c(-200, 400),
    xlim = c(0,365),
     xlab = "Date",
     ylab = "Water deficit in P-PET in mm",
     type = "l")
axis(1,at = axisticks, labels = months)

polygon(c(days, rev(days)), 
          c(clim.def.max.2014_2017, rev(clim.def.min.2014_2017)),
        col = "#AED6F1AA")
polygon(c(days, rev(days)), 
          c(clim.def.max.2018_2021, rev(clim.def.min.2018_2021)),
        col = "#F5B7B1AA")

Same plot as above, shading the between maximum and minimum of the periods: Blue is 2014-2017, red is 2018-2021.

Epot.yearly.mat <- matrix(h.data$Epot, nrow=365, ncol=length(years)-1 )
P.yearly.mat <- matrix(h.data$RSK, nrow=365, ncol=length(years)-1 )

Epot.yearly <- apply(Epot.yearly.mat,2,sum)
P.yearly <- apply(P.yearly.mat,2,sum)

(mean(Epot.yearly))

[1] 511.5383

(mean(P.yearly))

[1] 529.05

plot(h.years, P.yearly,
     pch = 19,
     xlab = "Hydrological year 2005/06 - 2020/21",
     ylim = c(0,900),
     ylab = "Flux in mm / hydrological year")
points(h.years, Epot.yearly,
       pch = 19,
       col = "red")

legend(x="topright",
       legend = c("precip", "Epot"),
       pch = c(1, 1),
       col = c("black", "red"))

       

lm.P <- lm(P.yearly ~ h.years)
(summary(lm.P))

Call:
lm(formula = P.yearly ~ h.years)

Residuals:
     Min       1Q   Median       3Q      Max 
-203.676  -96.138    7.434   80.568  200.481 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 19167.546  14054.444   1.364    0.194
h.years        -9.257      6.980  -1.326    0.206

Residual standard error: 128.7 on 14 degrees of freedom
Multiple R-squared:  0.1116,    Adjusted R-squared:  0.04815 
F-statistic: 1.759 on 1 and 14 DF,  p-value: 0.206

lm.Epot <- lm(Epot.yearly ~ h.years)
(summary(lm.Epot))

Call:
lm(formula = Epot.yearly ~ h.years)

Residuals:
   Min     1Q Median     3Q    Max 
-36.20 -16.78  -0.32   9.87  49.36 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -3259.587   2548.230  -1.279    0.222
h.years         1.873      1.266   1.480    0.161

Residual standard error: 23.34 on 14 degrees of freedom
Multiple R-squared:  0.1353,    Adjusted R-squared:  0.07351 
F-statistic:  2.19 on 1 and 14 DF,  p-value: 0.161

Annual PET significantly increases by 1.6 mm / year (due both the average temperature increase, and increase in daily temperature range). The average was 662 mm. Precipitation is more variable throughout the period 1992-2022. The average was 672 mm, no significant changes.

lm.Epot.yearly <- lm(Epot.yearly ~ h.years)
summary(lm.Epot.yearly)

Call:
lm(formula = Epot.yearly ~ h.years)

Residuals:
   Min     1Q Median     3Q    Max 
-36.20 -16.78  -0.32   9.87  49.36 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -3259.587   2548.230  -1.279    0.222
h.years         1.873      1.266   1.480    0.161

Residual standard error: 23.34 on 14 degrees of freedom
Multiple R-squared:  0.1353,    Adjusted R-squared:  0.07351 
F-statistic:  2.19 on 1 and 14 DF,  p-value: 0.161

lm.P.yearly <- lm(P.yearly ~ h.years)
summary(lm.P.yearly)

Call:
lm(formula = P.yearly ~ h.years)

Residuals:
     Min       1Q   Median       3Q      Max 
-203.676  -96.138    7.434   80.568  200.481 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 19167.546  14054.444   1.364    0.194
h.years        -9.257      6.980  -1.326    0.206

Residual standard error: 128.7 on 14 degrees of freedom
Multiple R-squared:  0.1116,    Adjusted R-squared:  0.04815 
F-statistic: 1.759 on 1 and 14 DF,  p-value: 0.206

def.yearly <- P.yearly - Epot.yearly
plot(h.years, def.yearly,
    xlab = "Hydrological year 1991/2 - 2021/22",
    ylab = "Water deficit or surplus in mm / hydrological year",
    main = "No 6265 Wusteritz")

lm.def.yearly <- lm(def.yearly ~ h.years)
(summary(lm.def.yearly))

Call:
lm(formula = def.yearly ~ h.years)

Residuals:
    Min      1Q  Median      3Q     Max 
-222.75 -114.93   11.86  104.20  190.66 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 22427.13   15444.26   1.452    0.169
h.years       -11.13       7.67  -1.451    0.169

Residual standard error: 141.4 on 14 degrees of freedom
Multiple R-squared:  0.1307,    Adjusted R-squared:  0.06864 
F-statistic: 2.105 on 1 and 14 DF,  p-value: 0.1688

The difference of the two (which is the meteorological surplus when positive or water deficit when negative), decreases, but not significantly.

plot(P.yearly, Epot.yearly,
     xlab = "Annual precipitation in mm / hydrological year",
     ylab = "Annual PET in mm / hydrological year")

lm.Epot.P <- lm(P.yearly ~ Epot.yearly)
summary(lm.Epot.P)

Call:
lm(formula = P.yearly ~ Epot.yearly)

Residuals:
    Min      1Q  Median      3Q     Max 
-119.34  -91.03  -19.77   63.15  253.68 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 2046.284    624.343   3.278   0.0055 **
Epot.yearly   -2.966      1.219  -2.433   0.0290 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 114.5 on 14 degrees of freedom
Multiple R-squared:  0.2971,    Adjusted R-squared:  0.2469 
F-statistic: 5.918 on 1 and 14 DF,  p-value: 0.02899

There is a significant relation between potential evapotranspiration and precipitation (calculated over the hydrological year). The peak of the precipitation is in summer. Hence, when precipitation is high, probably cloudiness is high and temperatures lower. In contrast, when precipitation is low, raindays less, sunshine hours increased, also the potential evaporation is high. T

Multiannual cumulative water deficit

Multiannual water deficit, starting the water balance in Octber 2005 and NOT resetting each year. The negative numbers can be intepreted as the water storage in the subsurface that would be required to compensate for the lack of precipitation.

h.data$cumdef <- cumsum(h.data$RSK - h.data$Epot) # longterm cumulated water deficit

plot(h.data$date, h.data$cumdef, 
#     ylim = c(-200,250),
     xlab = "Date",
     ylab = "P - PET since Oct 1 1991 in mm",
     main = "No 6265 Wusteritz")

Here we are looking at the surplus building up (site is energy limited before 2015) than reducing again (site is water limited after 2015). It is back to neutral in 2020-2021.

h.data$cumEpot <- cumsum(h.data$Epot) # longterm cumulated water deficit


plot(h.data$date, cumsum(h.data$RSK), 
     ylim = c(0,25000),
     xlab = "Date",
     ylab = "Cumulared flux since Oct 1 2005 in mm")
points(h.data$date, cumsum(h.data$Epot), 
#     ylim = c(-200,250),
     xlab = "date",
     col = "red")

legend(x = "topleft",          # Position
       legend = c("precipitation", "Epot"),  # Legend texts
       pch = c(1, 1),
       col = c("black", "red")
       ) # point colors

NA
NA

Precipitation analysis

QQplot.

q.p_daily.all <- quantile(data$RSK, seq(0.025,1,0.025))
q.p_daily.2014_2017 <- quantile(data_2014_2017$RSK, seq(0.025,1,0.025))
q.p_daily.2018_2021 <- quantile(data_2018_2021$RSK, seq(0.025,1,0.025))

q.p_daily.2020 <- quantile(data_2020$RSK, seq(0.025,1,0.025))
q.p_daily.2021 <- quantile(data_2021$RSK, seq(0.025,1,0.025))


plot(q.p_daily.all, q.p_daily.2014_2017, 
     xlab = "Quantiles of daily precipitation 2005-2021",
     ylab = "Quantiles of daily precipitation in the observation period",
     xlim = c(0,100),
     ylim = c(0,100))
points(q.p_daily.all, q.p_daily.2018_2021, 
        col = "red")

abline(a=0, b=1, ':')

NAs introduced by coercion

plot(q.p_daily.2014_2017, q.p_daily.2018_2021, 
     ylim = c(0,100),
     xlab = "Quantiles of daily precipitation 2014-2017",
     ylab = "Quantiles of daily precipitation 2018-2021")
abline(a=0, b=1, ':')

NAs introduced by coercion

plot(q.p_daily.2020, q.p_daily.2021, 
     ylim = c(0,100),
     xlab = "Quantiles of daily precipitation 2020",
     ylab = "Quantiles of daily precipitation 2021")
abline(a=0, b=1, ':')

NAs introduced by coercion

plot(q.p_daily.2020, q.p_daily.2021, 
     ylim = c(0,100),
     xlab = "Quantiles of daily precipitation 2020",
     ylab = "Quantiles of daily precipitation 2021")
abline(a=0, b=1, ':')

NAs introduced by coercion

Table_QP <- as.data.frame(cbind(seq(0.025,1,0.025), q.p_daily.2020, q.p_daily.2021))
names(Table_QP) <- c("P", "QP_2020", "QP_2021")

write.csv(Table_QP, "Table_QP_6265.csv")

Number of days with given rainfall

Check whether certain types of daily precipitation were more or less common along the time series. For this I use the entire time series, to be able to see trends (if applicable).

I will test for

• number of days per year without precipitation (threshold 0.1 mm)

• days surpassing small precipitation (threshold set to 5 mm)

Looking for extreme rainfall requires a dataset in hourly or six hourly resolution, and cannot be done using this daily dataset. But the data is available from DWD.

# intialize
ndays = array(data=NA,dim = nyears, dimnames = NULL)
noraindays = array(data=NA,dim = nyears, dimnames = NULL)
strongraindays = array(data=NA,dim = nyears, dimnames = NULL)
P = array(data=NA,dim = nyears, dimnames = NULL)
P_strong = array(data=NA,dim = nyears, dimnames = NULL)

P_thresh = 5

#i=1
for (i in 1:nyears) {
data_yr = subset(data, format(data$date,'%Y') == as.character(years[i]))
noraindays[i] = length(which(data_yr$RSK<0.1))
raindays = nrow(data_yr)-noraindays
strongraindays[i] = length(which(data_yr$RSK>P_thresh ))
P[i] <- sum(data_yr$RSK)
P_strong[i] <- sum(data_yr$RSK[which(data_yr$RSK>P_thresh )])
ndays[i] = nrow(data_yr)

}

Annual precipitation time series

par(mfrow=c(1,3))
plot(years, P,
     ylab = "annual precpitation in mm/a",
     main = "all days")

plot(years, P-P_strong,
     ylab = "annual precpitation < 5 mm in mm/a",
     main = "small precipitation")

plot(years, P_strong,
     ylab = "annual precpitation > 5mm in mm/a",
     main = "precipitation > 5 mm/d")

NA
NA

lm.p_yr <-lm(P~years) # all years
summary(lm.p_yr)

Call:
lm(formula = P ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-210.52  -71.01  -23.64   94.67  214.44 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 15950.753  11682.324   1.365    0.192
years          -7.661      5.803  -1.320    0.207

Residual standard error: 117.2 on 15 degrees of freedom
Multiple R-squared:  0.1041,    Adjusted R-squared:  0.04435 
F-statistic: 1.743 on 1 and 15 DF,  p-value: 0.2066

yrs_2014_2021 = c(10:17)

lm.p_yr.2014_2021 <-lm(P[yrs_2014_2021]~years[yrs_2014_2021]) # 2014-2021 only
summary(lm.p_yr.2014_2021)

Call:
lm(formula = P[yrs_2014_2021] ~ years[yrs_2014_2021])

Residuals:
     Min       1Q   Median       3Q      Max 
-175.906  -63.965   -5.087   90.648  133.031 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)
(Intercept)          26230.63   34872.00   0.752    0.480
years[yrs_2014_2021]   -12.76      17.28  -0.738    0.488

Residual standard error: 112 on 6 degrees of freedom
Multiple R-squared:  0.0833,    Adjusted R-squared:  -0.06948 
F-statistic: 0.5452 on 1 and 6 DF,  p-value: 0.4881

No obvious temporal trend in annual precipitation.

P_weak <- P-P_strong
lm.p_weak_yr <-lm(P_weak~years) # all years
summary(lm.p_weak_yr)

Call:
lm(formula = P_weak ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-27.603  -7.441   0.176  11.206  32.071 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 4071.7206  1664.0533   2.447   0.0272 *
years         -1.9265     0.8267  -2.330   0.0342 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 16.7 on 15 degrees of freedom
Multiple R-squared:  0.2658,    Adjusted R-squared:  0.2169 
F-statistic: 5.431 on 1 and 15 DF,  p-value: 0.03415

yrs_2014_2018 = c(10:17)

lm.pw_eak_yr.2014_2021 <-lm(P_weak[yrs_2014_2021]~years[yrs_2014_2021]) # 2014-2021 only
summary(lm.pw_eak_yr.2014_2021)

Call:
lm(formula = P_weak[yrs_2014_2021] ~ years[yrs_2014_2021])

Residuals:
    Min      1Q  Median      3Q     Max 
-29.114 -12.377  -2.025  13.409  29.164 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)
(Intercept)          6888.257   6911.342   0.997    0.357
years[yrs_2014_2021]   -3.321      3.426  -0.970    0.370

Residual standard error: 22.2 on 6 degrees of freedom
Multiple R-squared:  0.1355,    Adjusted R-squared:  -0.008638 
F-statistic: 0.9401 on 1 and 6 DF,  p-value: 0.3697

Decrease of low intensity precipitation (P<5 mm/d) part of annual precipitation over the entire period, but not significant over the shorter period.

lm.p_strong_yr <-lm(P_strong~years) # all years
summary(lm.p_strong_yr)

Call:
lm(formula = P_strong ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-203.50  -85.06  -16.20   88.14  203.24 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 11879.032  11286.865   1.052    0.309
years          -5.734      5.607  -1.023    0.323

Residual standard error: 113.3 on 15 degrees of freedom
Multiple R-squared:  0.06518,   Adjusted R-squared:  0.002863 
F-statistic: 1.046 on 1 and 15 DF,  p-value: 0.3227

summary(lm(P_strong/P~years))

Call:
lm(formula = P_strong/P ~ years)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.16566 -0.06364  0.01439  0.05841  0.09892 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.589820   8.702436   0.413    0.686
years       -0.001476   0.004323  -0.341    0.737

Residual standard error: 0.08732 on 15 degrees of freedom
Multiple R-squared:  0.007714,  Adjusted R-squared:  -0.05844 
F-statistic: 0.1166 on 1 and 15 DF,  p-value: 0.7375

yrs_2014_2018 = c(10:17)

lm.p_strong_yr.2014_2021 <-lm(P_strong[yrs_2014_2021]~years[yrs_2014_2021]) # 2014-2021 only
summary(lm.p_strong_yr.2014_2021)

Call:
lm(formula = P_strong[yrs_2014_2021] ~ years[yrs_2014_2021])

Residuals:
     Min       1Q   Median       3Q      Max 
-146.792  -77.819    9.988   90.235  103.867 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)          19342.375  32578.290   0.594    0.574
years[yrs_2014_2021]    -9.442     16.148  -0.585    0.580

Residual standard error: 104.6 on 6 degrees of freedom
Multiple R-squared:  0.05391,   Adjusted R-squared:  -0.1038 
F-statistic: 0.3419 on 1 and 6 DF,  p-value: 0.5801

plot(years, P_strong/P)

plot(years, strongraindays/raindays)

Days per year without precipitation.

plot(years, noraindays)

lm.norain <- lm(noraindays~years) # all years
summary(lm.norain)

Call:
lm(formula = noraindays ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-23.137 -10.770  -7.892   7.966  30.740 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -2062.0441  1521.3875  -1.355    0.195
years           1.1225     0.7558   1.485    0.158

Residual standard error: 15.27 on 15 degrees of freedom
Multiple R-squared:  0.1282,    Adjusted R-squared:  0.0701 
F-statistic: 2.206 on 1 and 15 DF,  p-value: 0.1582

lm.norain.2014_2021 <- lm(noraindays[yrs_2014_2018]~years[yrs_2014_2018]) # all years
summary(lm.norain.2014_2021)

Call:
lm(formula = noraindays[yrs_2014_2018] ~ years[yrs_2014_2018])

Residuals:
     Min       1Q   Median       3Q      Max 
-21.1786  -7.6518  -0.4464   4.5625  29.9286 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)          -7651.714   5706.911  -1.341    0.229
years[yrs_2014_2018]     3.893      2.829   1.376    0.218

Residual standard error: 18.33 on 6 degrees of freedom
Multiple R-squared:  0.2399,    Adjusted R-squared:  0.1132 
F-statistic: 1.894 on 1 and 6 DF,  p-value: 0.2179

Days with some rain

weakraindays = raindays - strongraindays
plot(years, weakraindays)

lm.weakrain <- lm(weakraindays~years)
summary(lm.weakrain)

Call:
lm(formula = weakraindays ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.005  -8.936   4.034   7.520  16.015 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 1134.1029  1279.9382   0.886    0.390
years         -0.4951     0.6358  -0.779    0.448

Residual standard error: 12.84 on 15 degrees of freedom
Multiple R-squared:  0.03885,   Adjusted R-squared:  -0.02523 
F-statistic: 0.6063 on 1 and 15 DF,  p-value: 0.4483

Days with substantial rain

plot(years, strongraindays)

lm.strongrain <- lm(strongraindays~years)
summary(lm.strongrain)

Call:
lm(formula = strongraindays ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.274  -6.118   1.137   5.608   8.627 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1292.9412   650.5286   1.988   0.0654 .
years         -0.6275     0.3232  -1.942   0.0712 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 6.528 on 15 degrees of freedom
Multiple R-squared:  0.2008,    Adjusted R-squared:  0.1476 
F-statistic:  3.77 on 1 and 15 DF,  p-value: 0.07121

table_PET_P <- as.data.frame(cbind(h.years,  P.yearly, Epot.yearly))
colnames(table_PET_P) <- c("years", "P", "PET")
write.csv(table_PET_P, file="P_PET_16yrs_6265.csv")

Ptable <- as.data.frame(cbind(years, P, P-P_strong, P_strong, noraindays, weakraindays, strongraindays))

colnames(Ptable) <- c("year", "P","P_small", "P_subst", "ndays_norain","ndays_weakrain","ndays_substrain")
write.csv(Ptable, file <- "Ptable_16yrs_6265.csv")