Hainich 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 6305 Mühlhausen/Thüringen-Görmar (01.12.2004 - current) Kläranlage

filename_mhlsn = "tageswerte_KL_06305_20041201_20221231_hist/produkt_klima_tag_20041201_20221231_06305.txt" # Mühlhausen - nearest dwd precip station

There were 20 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] 31 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). 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

nyears = length(years)

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, no missing values

integer(0)

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

integer(0)

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

integer(0)

#(a <- which(data$SDK == -999)) # daily sunshine duration, no missing values
#a <- which(data$NM == -999) # daily cloudcover, very many missing values, dó not use this

(b <- which(data$RSK == -999)) # daily precipitation duration, one missing value. Replace with zero

[1] 5671

data$date[b]

[1] "2020-07-15"

data$RSK[b] <- 0

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)-1]

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.

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

# 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),] 

# 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.4079412  0.4441176  2.0923529 

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

       10%        50%        90% 
-1.2870588  0.5073529  2.4188235 

# temperature anomalies
data$Tclim <- NA
ind <- which(data$year>2004)
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.5431  -2.2387   0.0319   2.3763  11.4615 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -1.100e+02  1.855e+01  -5.929 3.22e-09 ***
as.numeric(year)  5.465e-02  9.217e-03   5.929 3.21e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.557 on 6203 degrees of freedom
Multiple R-squared:  0.005635,  Adjusted R-squared:  0.005475 
F-statistic: 35.15 on 1 and 6203 DF,  p-value: 3.215e-09

Both periods were warmer than the entire reference periods from 2005-2022, 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.06 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")

bins = seq(1,10,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)

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("clim", "2014-2017", "2018-2021"),  # Legend texts
       pch = c(1, 19, 19), # point types
       col = c("black", "blue", "red"), # point colors
       )           

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

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

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_6305.csv")

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,100),
     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 )

data length [5475] is not a sub-multiple or multiple of the number of columns [16]

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 <- which(h.years>=2014 & h.years<=2017)

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 
-82.504 -32.596  -4.361  20.374 143.675 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)
(Intercept)            -8581.370   6188.733  -1.387    0.187
years[2:length(years)]     4.372      3.074   1.422    0.177

Residual standard error: 56.67 on 14 degrees of freedom
Multiple R-squared:  0.1263,    Adjusted R-squared:  0.06385 
F-statistic: 2.023 on 1 and 14 DF,  p-value: 0.1768

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 
-28.995 -12.252  -5.939   9.934  51.333 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)
(Intercept)            -2051.748   2272.299  -0.903    0.382
years[2:length(years)]     1.049      1.129   0.930    0.368

Residual standard error: 20.81 on 14 degrees of freedom
Multiple R-squared:  0.05817,   Adjusted R-squared:  -0.009108 
F-statistic: 0.8646 on 1 and 14 DF,  p-value: 0.3682

#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 increases in the recent years, and arguably the variation in the water balance as well. The increase is not significant however.

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(-100, 150),
     xlab = "Date",
     ylab = "Water deficit in 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 negative balance in summer, whhn ET is larger. The most striking difference between the first and second period is that the surplus in winter is smaller in the second period. The deficit does not necessarily develop earlier in summer though (the red curve is still positive in July, while the blue one is already negative at this point).

#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, 200),
    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.

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

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

data length [5475] is not a sub-multiple or multiple of the number of columns [16]

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

data length [5475] is not a sub-multiple or multiple of the number of columns [16]

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

(mean(Epot.yearly))

[1] 547.6599

(mean(P.yearly))

[1] 502.55

plot(h.years, P.yearly,
     pch = 19,
     xlab = "Hydrological year 2005/06 - 2020/21",
     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 
-96.74 -71.16 -22.68  76.34 135.55 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 8772.113   8730.571   1.005    0.332
h.years       -4.107      4.336  -0.947    0.360

Residual standard error: 79.95 on 14 degrees of freedom
Multiple R-squared:  0.06023,   Adjusted R-squared:  -0.006902 
F-statistic: 0.8972 on 1 and 14 DF,  p-value: 0.3596

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-34.560 -11.125  -0.157  14.449  45.870 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -3965.249   2264.320  -1.751   0.1018  
h.years         2.241      1.125   1.993   0.0661 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 20.74 on 14 degrees of freedom
Multiple R-squared:  0.221, Adjusted R-squared:  0.1654 
F-statistic: 3.972 on 1 and 14 DF,  p-value: 0.06612

Annual PET slightly increases (due both the average temperature increase, and increase in daily temperature range) not significantly, but with the trend. The average was 548 mm. Precipitation is more variable throughout the period 2005/6-2021/22. The average was 501 mm there was no trend.

def.yearly <- P.yearly - Epot.yearly

Error: object 'P.yearly' not found

Also the difference of the two (which is the water deficit) does not give a significant relation. Overall, since the hydrological year 2017/18 the water deficit has not been positive, therefore adding up to a substantial water deficit at the end of 2021 seen below.

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 
-96.14 -46.27 -21.75  43.70 155.65 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) 1583.6581   425.2127   3.724  0.00227 **
Epot.yearly   -1.9741     0.7758  -2.545  0.02336 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 68.2 on 14 degrees of freedom
Multiple R-squared:  0.3162,    Adjusted R-squared:  0.2674 
F-statistic: 6.475 on 1 and 14 DF,  p-value: 0.02336

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. The combination leads to a substantial build up of water deficit that cannot be cured with one good year of precipitation.

Multiannual cumulative water deficit

Most impressive is the 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 2005 in mm",
     main = "No 6305 Muehlhausen")

Here we are looking at the deficit building up over several years. Note this was calculated using the potential evapotranspiration, not the real one. Regardless, the meteolorogical water deficit defenitiely built up, specifically since 2014. While the water deficit was rather balanced until the year 2014, the meteorological water balance went strongly negative since 2015, with some recovery in spring 2018 and an even stronger drop in summer 2018. The deficit at the end of 2021 corresponds roughly to one year of precipitation.

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


plot(h.data$date, cumsum(h.data$RSK), 
     ylim = c(0,9500),
     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

The precipitation and potential evaporation were almost equal until roughly 2015, when they started to deviate for several years in a row. The plot shows that the discrepancy built up since 2015, but was much more pronounced after 2018.

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))

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(log(q.p_daily.2014_2017+1), log(q.p_daily.2018_2021+1), 
     ylim = c(0,3),
     xlim = c(0,3),
     xlab = "Quantiles of daily precipitation 2014-2017",
     ylab = "Quantiles of daily precipitation 2018-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.2014_2017, q.p_daily.2018_2021))
names(Table_QP) <- c("P", "QP_2014_2017", "QP_2018_2021")

write.csv(Table_QP, "Table_QP_6305.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 
-139.367  -65.589   -7.811   50.507  155.185 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 7431.272   8308.957   0.894    0.385
years         -3.441      4.128  -0.834    0.418

Residual standard error: 83.37 on 15 degrees of freedom
Multiple R-squared:  0.04427,   Adjusted R-squared:  -0.01944 
F-statistic: 0.6948 on 1 and 15 DF,  p-value: 0.4176

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 
-133.244  -27.769    5.306   32.660  122.169 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)
(Intercept)           -755.3821 24724.7872  -0.031    0.977
years[yrs_2014_2021]     0.6131    12.2552   0.050    0.962

Residual standard error: 79.42 on 6 degrees of freedom
Multiple R-squared:  0.000417,  Adjusted R-squared:  -0.1662 
F-statistic: 0.002503 on 1 and 6 DF,  p-value: 0.9617

No obvious temporal trend in annual precipitation.

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

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

Residuals:
      Min        1Q    Median        3Q       Max 
-0.092882 -0.047994 -0.002891  0.039094  0.158030 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  7.420787   6.861925   1.081    0.297
years       -0.003494   0.003409  -1.025    0.322

Residual standard error: 0.06885 on 15 degrees of freedom
Multiple R-squared:  0.06546,   Adjusted R-squared:  0.00316 
F-statistic: 1.051 on 1 and 15 DF,  p-value: 0.3216

lm.weak_raindays <-lm(weakraindays~years) # all years
summary(lm.weak_raindays)

Call:
lm(formula = weakraindays ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.029  -6.824  -3.618   8.147  17.559 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 2565.5000  1057.7118   2.426   0.0284 *
years         -1.2059     0.5254  -2.295   0.0366 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.61 on 15 degrees of freedom
Multiple R-squared:  0.2599,    Adjusted R-squared:  0.2105 
F-statistic: 5.267 on 1 and 15 DF,  p-value: 0.03658

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 
-33.168 -15.161  -8.182  17.986  35.904 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)
(Intercept)          4647.496   8238.362   0.564    0.593
years[yrs_2014_2021]   -2.214      4.083  -0.542    0.607

Residual standard error: 26.46 on 6 degrees of freedom
Multiple R-squared:  0.04672,   Adjusted R-squared:  -0.1122 
F-statistic: 0.294 on 1 and 6 DF,  p-value: 0.6072

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 
-127.52  -54.25   18.68   28.69  148.96 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 1770.0581  7998.7252   0.221    0.828
years         -0.7238     3.9735  -0.182    0.858

Residual standard error: 80.26 on 15 degrees of freedom
Multiple R-squared:  0.002207,  Adjusted R-squared:  -0.06431 
F-statistic: 0.03318 on 1 and 15 DF,  p-value: 0.8579

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 
-100.076  -34.065    3.583   30.440   90.651 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)          -5402.879  19812.165  -0.273    0.794
years[yrs_2014_2021]     2.827      9.820   0.288    0.783

Residual standard error: 63.64 on 6 degrees of freedom
Multiple R-squared:  0.01363,   Adjusted R-squared:  -0.1508 
F-statistic: 0.0829 on 1 and 6 DF,  p-value: 0.7831

plot(years, P_strong/P)

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

plot(years, strongraindays/raindays)

summary(lm(P_strong/P~years))

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

Residuals:
      Min        1Q    Median        3Q       Max 
-0.158030 -0.039094  0.002891  0.047994  0.092882 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.420787   6.861925  -0.936    0.364
years        0.003494   0.003409   1.025    0.322

Residual standard error: 0.06885 on 15 degrees of freedom
Multiple R-squared:  0.06546,   Adjusted R-squared:  0.00316 
F-statistic: 1.051 on 1 and 15 DF,  p-value: 0.3216

summary(lm(strongraindays/raindays~years))

Call:
lm(formula = strongraindays/raindays ~ years)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.058075 -0.026681  0.001203  0.018081  0.061480 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.547222   3.582053  -0.711    0.488
years        0.001356   0.001779   0.762    0.458

Residual standard error: 0.03594 on 15 degrees of freedom
Multiple R-squared:  0.03728,   Adjusted R-squared:  -0.0269 
F-statistic: 0.5809 on 1 and 15 DF,  p-value: 0.4578

Days per year without precipitation.

summary(lm.norain)

Call:
lm(formula = noraindays ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-16.360 -10.135  -2.973   5.703  26.091 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -2206.8309  1281.8557  -1.722   0.1057  
years           1.1936     0.6368   1.874   0.0805 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.86 on 15 degrees of freedom
Multiple R-squared:  0.1898,    Adjusted R-squared:  0.1358 
F-statistic: 3.514 on 1 and 15 DF,  p-value: 0.08048

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 
-16.029  -6.824  -3.618   8.147  17.559 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 2565.5000  1057.7118   2.426   0.0284 *
years         -1.2059     0.5254  -2.295   0.0366 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.61 on 15 degrees of freedom
Multiple R-squared:  0.2599,    Adjusted R-squared:  0.2105 
F-statistic: 5.267 on 1 and 15 DF,  p-value: 0.03658

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 
-10.061  -4.902  -1.073   4.951  12.049 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)   6.33088  715.85100   0.009    0.993
years         0.01225    0.35561   0.034    0.973

Residual standard error: 7.183 on 15 degrees of freedom
Multiple R-squared:  7.917e-05, Adjusted R-squared:  -0.06658 
F-statistic: 0.001188 on 1 and 15 DF,  p-value: 0.973

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_6305.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_6305.csv")