SESO 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 3289 Schmierits-Weltwitz (01.10.1972 - current) near SESO - Data missing - start analysis ti fit Hai, which only starts in 2005

filename = "tageswerte_KL_03289_19721001_20221231_hist/produkt_klima_tag_19721001_20221231_03289.txt" # Schmierits-Weltwitz

There were 28 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:49] 92 365 365 365 366 365 365 365 366 365 ...
  values : chr [1:49] "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979" "1980" "1981" "1982" "1983" "1984" ...

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

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

days_per_year[[1]]

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

Selected years are complete.

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

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

[1] 5334

(data$date[b]) # need to check surrounging stations, e.g. Stadtroda For now take as zero

[1] "2019-08-09"

data$RSK[b] <- 0

Data for temperature is complete. Missing precipitation in August 2019. Need to double check. For now set to zero.

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 and other variables changed between

• 2014-2017 and

• 2018-2021.

Since the time series of SESO is a bit shorter (2016-2021), create an alterative subset, to be able to look only at some of the analysis for this smaller range.

• 2016-2018 and

• 2019-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),] 

data_2016_2018 <- data[which(data$year>=2016 & data$year<=2018),] 
data_2019_2021 <- data[which(data$year>=2019 & 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.7426471  0.4970588  2.4404412 

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

      10%       50%       90% 
-1.336029  0.425000  2.625735 

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

number of items to replace is not a multiple of replacement length

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 
-18.5489  -2.7984   0.0051   2.7944  12.3997 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -117.16155   23.14910  -5.061 4.29e-07 ***
as.numeric(year)    0.05821    0.01150   5.063 4.26e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.051 on 5842 degrees of freedom
  (365 observations deleted due to missingness)
Multiple R-squared:  0.004368,  Adjusted R-squared:  0.004198 
F-statistic: 25.63 on 1 and 5842 DF,  p-value: 4.261e-07

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 2014-2017
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 2018-2021
DeltaT.clim.2018_2021 <- aggregate(DeltaT ~ day+month, data = data_2018_2021, FUN = mean)
points(DeltaT.clim.2018_2021$DeltaT, col="red")

# averge over periods 2016-2018
DeltaT.clim.2016_2018 <- aggregate(DeltaT ~ day+month, data = data_2016_2018, FUN = mean)
#points(DeltaT.clim.2016_2018$DeltaT, col="blue")

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

bins = seq(0.5,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)
h.DeltaT_daily.2016_2018 <- hist(DeltaT.clim.2016_2018$DeltaT, breaks = bins, plot = FALSE)
h.DeltaT_daily.2019_2021 <- hist(DeltaT.clim.2019_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
       )           

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.2016_2018$mids, cumsum(h.DeltaT_daily.2016_2018$density)/4, 
       pch = 19, col = "blue")

points(h.DeltaT_daily.2019_2021$mids, cumsum(h.DeltaT_daily.2019_2021$density)/4, 
       pch = 19, col = "red")
legend(x = "bottomright",          # Position
       legend = c("clim", "2016-2018", "2019-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_3289.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)

package 'tidyr' was built under R version 4.0.5

# 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 
-95.92 -35.84 -12.64  51.49 137.85 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)
(Intercept)            2293.103   7661.781   0.299    0.769
years[2:length(years)]   -1.033      3.805  -0.271    0.790

Residual standard error: 70.16 on 14 degrees of freedom
Multiple R-squared:  0.005234,  Adjusted R-squared:  -0.06582 
F-statistic: 0.07366 on 1 and 14 DF,  p-value: 0.79

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 
-25.337 -15.246  -6.222  11.149  43.015 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)
(Intercept)             265.6458  2381.2035   0.112    0.913
years[2:length(years)]   -0.1044     1.1826  -0.088    0.931

Residual standard error: 21.81 on 14 degrees of freedom
Multiple R-squared:  0.0005562, Adjusted R-squared:  -0.07083 
F-statistic: 0.007791 on 1 and 14 DF,  p-value: 0.9309

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

(mean(P.yearly))

[1] 665.9937

plot(h.years, P.yearly,
     pch = 19,
     xlab = "Hydrological year 2015/16 - 2021/22",
     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 
-184.190  -72.007   -7.682   73.639  191.949 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 15023.729  12721.873   1.181    0.257
h.years        -7.131      6.318  -1.129    0.278

Residual standard error: 116.5 on 14 degrees of freedom
Multiple R-squared:  0.08339,   Adjusted R-squared:  0.01792 
F-statistic: 1.274 on 1 and 14 DF,  p-value: 0.278

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

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

Residuals:
    Min      1Q  Median      3Q     Max 
-55.072  -9.141   4.969  22.300  35.472 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) -5117.087   3197.800  -1.600   0.1319  
h.years         2.877      1.588   1.811   0.0916 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 29.28 on 14 degrees of freedom
Multiple R-squared:  0.1899,    Adjusted R-squared:  0.132 
F-statistic: 3.281 on 1 and 14 DF,  p-value: 0.0916

table_3289 <- c(h.years,  P.yearly, Epot.yearly)
write.csv(table_3289, file="P_PET_30yrs_3289.csv")

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.

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 3289 Schmieritz-Weltwitz")

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 
-207.856  -86.016   -9.599   66.005  236.317 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 20140.817  14693.334   1.371    0.192
h.years       -10.007      7.297  -1.371    0.192

Residual standard error: 134.6 on 14 degrees of freedom
Multiple R-squared:  0.1184,    Adjusted R-squared:  0.05545 
F-statistic: 1.881 on 1 and 14 DF,  p-value: 0.1918

Also the difference of the two (which is the meteorological surplus when positive or water deficit when negative) also decreases significantly by 5.6 mm/year.

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 
-163.578  -53.982   -0.215   24.785  237.514 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) 2155.4621   545.3838   3.952  0.00145 **
Epot.yearly   -2.2065     0.8071  -2.734  0.01615 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 98.26 on 14 degrees of freedom
Multiple R-squared:  0.348, Adjusted R-squared:  0.3015 
F-statistic: 7.474 on 1 and 14 DF,  p-value: 0.01615

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 3289 Schmieritz-Weltwitz")

Here we are looking at the surplus building up over several years than disappearing again.

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

QQplot.

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

There were 50 or more warnings (use warnings() to see the first 50)

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 
-179.55 -104.16  -36.83   66.32  235.59 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 12680.607  11987.351   1.058    0.307
years          -5.970      5.955  -1.003    0.332

Residual standard error: 120.3 on 15 degrees of freedom
Multiple R-squared:  0.0628,    Adjusted R-squared:  0.0003183 
F-statistic: 1.005 on 1 and 15 DF,  p-value: 0.332

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 
-173.79  -23.39   21.09   54.62   78.09 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)          -8199.161  28258.937  -0.290    0.781
years[yrs_2014_2021]     4.374     14.007   0.312    0.765

Residual standard error: 90.78 on 6 degrees of freedom
Multiple R-squared:  0.01599,   Adjusted R-squared:  -0.148 
F-statistic: 0.09751 on 1 and 6 DF,  p-value: 0.7654

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 
-44.821  -8.814  -2.297  20.382  27.672 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 3421.610   2170.207   1.577    0.136
years         -1.597      1.078  -1.481    0.159

Residual standard error: 21.78 on 15 degrees of freedom
Multiple R-squared:  0.1276,    Adjusted R-squared:  0.0694 
F-statistic: 2.193 on 1 and 15 DF,  p-value: 0.1593

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 
-37.499 -15.698  -0.501  17.806  36.596 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)          -4538.796   8754.746  -0.518    0.623
years[yrs_2014_2021]     2.348      4.339   0.541    0.608

Residual standard error: 28.12 on 6 degrees of freedom
Multiple R-squared:  0.04651,   Adjusted R-squared:  -0.1124 
F-statistic: 0.2927 on 1 and 6 DF,  p-value: 0.608

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 
-141.015  -92.271    7.991   51.532  234.097 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  9258.997  11158.878   0.830    0.420
years          -4.374      5.543  -0.789    0.442

Residual standard error: 112 on 15 degrees of freedom
Multiple R-squared:  0.03984,   Adjusted R-squared:  -0.02417 
F-statistic: 0.6225 on 1 and 15 DF,  p-value: 0.4424

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 
-136.29  -48.55   38.50   45.45   63.13 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)          -3660.364  23756.482  -0.154    0.883
years[yrs_2014_2021]     2.026     11.775   0.172    0.869

Residual standard error: 76.31 on 6 degrees of freedom
Multiple R-squared:  0.004911,  Adjusted R-squared:  -0.1609 
F-statistic: 0.02961 on 1 and 6 DF,  p-value: 0.869

summary(lm(P_weak/P~ years))

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

Residuals:
     Min       1Q   Median       3Q      Max 
-0.08636 -0.01678 -0.00061  0.03590  0.06078 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.6461596  4.6297328   0.140    0.891
years       -0.0001623  0.0022999  -0.071    0.945

Residual standard error: 0.04646 on 15 degrees of freedom
Multiple R-squared:  0.0003319, Adjusted R-squared:  -0.06631 
F-statistic: 0.004981 on 1 and 15 DF,  p-value: 0.9447

Days per year without precipitation.

summary(lm.norain)

Call:
lm(formula = noraindays ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-22.235  -5.691  -1.029   3.779  26.368 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -1559.2794  1319.2090  -1.182    0.256
years           0.8676     0.6553   1.324    0.205

Residual standard error: 13.24 on 15 degrees of freedom
Multiple R-squared:  0.1046,    Adjusted R-squared:  0.04494 
F-statistic: 1.753 on 1 and 15 DF,  p-value: 0.2053

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 
-19.743  -7.404   1.397   4.228  18.919 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  479.0221  1069.1685   0.448    0.661
years         -0.1691     0.5311  -0.318    0.755

Residual standard error: 10.73 on 15 degrees of freedom
Multiple R-squared:  0.006714,  Adjusted R-squared:  -0.05951 
F-statistic: 0.1014 on 1 and 15 DF,  p-value: 0.7546

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.0221  -6.8162   0.2794   6.0882  16.7868 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1445.2574   788.3681   1.833   0.0867 .
years         -0.6985     0.3916  -1.784   0.0947 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.911 on 15 degrees of freedom
Multiple R-squared:  0.175, Adjusted R-squared:   0.12 
F-statistic: 3.181 on 1 and 15 DF,  p-value: 0.09473

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