Precipitation analysis

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

Chose station No 6305 Mühlhausen/Thüringen-Görmar

filename_mhlsn = "tageswerte_RR_06305_20041201_20221231_hist/produkt_nieder_tag_20041201_20221231_06305.txt" # Mühlhausen - nearest dwd precip station

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

data=read.csv(filename_mhlsn, header = TRUE, sep=";" )

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 <- as.numeric(format(data$date,'%Y'))

(days_per_year <- rle(data$year))

Run Length Encoding
  lengths: int [1:19] 31 365 365 365 366 365 365 365 366 365 ...
  values : num [1:19] 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 ...

One year is incomplete, start with the first complete year, whcih is 2005

years <- years[-1]
nyears = length(years)

remove Feb 29

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

Basic features of the rainfall dataset

Check whether the distribution of rainfall changed between

• 2014-2017 and

• 2018-2021.

I chose two intervals of equal size, therefore not covereing the period before 2014.

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

maxP = max(data$RS)

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

bins <- c(0:ceiling(maxP))

h.p_daily.2014_2017 <- hist(data_2014_2017$RS, breaks = bins, plot = FALSE)
h.p_daily.2018_2021 <- hist(data_2018_2021$RS, breaks = bins, plot = FALSE)

plot(h.p_daily.2014_2017$mids, log(h.p_daily.2014_2017$density), 
     xlab = "Daily precipitation in mm/d",
     ylab = "log(occurence frequency")
points(h.p_daily.2018_2021$mids, log(h.p_daily.2018_2021$density), pch = 19)

legend(x = "topright",          # Position
       legend = c("2014-2017", "2018-2021"),  # Legend texts
       pch = c(1, 19))           # point types

It is hard to judge. There is a tendency of smaller events being smaller in 2018-2021 compared to the period 2014-2017. Another look at this allows the QQplot.

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

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

NAs introduced by coercion

And zoom in a bit on the small numbers by plotting log(quantile) instead of the quantiles.

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

NAs introduced by coercion

It looks like the small events were less common, but the larger increasingly more common in the later period (2018-2021) compared to the first period (2014-2017). It is a small shift.

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)
snowdays= array(data=NA,dim = nyears, dimnames = NULL)
P_snow <- 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$RS<0.1))
raindays = nrow(data_yr)-noraindays
strongraindays[i] = length(which(data_yr$RS>P_thresh ))
P[i] <- sum(data_yr$RS)
P_strong[i] <- sum(data_yr$RS[which(data_yr$RS>P_thresh )])
ndays[i] = nrow(data_yr)

snowdays[i] <-  length(which(data_yr$RSF==7))
P_snow[i] <- sum(data_yr$RS[which(data_yr$RSF==7 )])

}

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

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

Call:
lm(formula = P ~ years)

Residuals:
     Min       1Q   Median       3Q      Max 
-137.282  -59.227   -8.633   49.602  154.279 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 7978.097   7391.305   1.079    0.296
years         -3.713      3.671  -1.011    0.327

Residual standard error: 80.8 on 16 degrees of freedom
Multiple R-squared:  0.06009,   Adjusted R-squared:  0.001346 
F-statistic: 1.023 on 1 and 16 DF,  p-value: 0.3269

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~years) # all years
summary(lm.p_weak_yr)

Call:
lm(formula = P_weak ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-33.982 -15.370  -3.298  12.281  38.070 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept) 5295.065   2028.562   2.610   0.0189 *
years         -2.535      1.007  -2.516   0.0229 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 22.18 on 16 degrees of freedom
Multiple R-squared:  0.2835,    Adjusted R-squared:  0.2387 
F-statistic:  6.33 on 1 and 16 DF,  p-value: 0.02293

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 
-128.579  -50.189    2.456   32.223  147.444 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 2683.033   7128.398   0.376    0.712
years         -1.178      3.540  -0.333    0.744

Residual standard error: 77.93 on 16 degrees of freedom
Multiple R-squared:  0.006871,  Adjusted R-squared:  -0.0552 
F-statistic: 0.1107 on 1 and 16 DF,  p-value: 0.7437

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

No trend in strong precipitation.

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 
-16.4585  -9.8265  -0.5272   6.0900  24.8191 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -2543.7148  1156.9145  -2.199   0.0430 *
years           1.3612     0.5746   2.369   0.0308 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12.65 on 16 degrees of freedom
Multiple R-squared:  0.2597,    Adjusted R-squared:  0.2134 
F-statistic: 5.612 on 1 and 16 DF,  p-value: 0.03075

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 
-15.3333  -8.6786  -0.7143   3.0357  26.9524 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)          -6045.143   4394.505  -1.376    0.218
years[yrs_2014_2018]     3.095      2.178   1.421    0.205

Residual standard error: 14.12 on 6 degrees of freedom
Multiple R-squared:  0.2518,    Adjusted R-squared:  0.1271 
F-statistic: 2.019 on 1 and 6 DF,  p-value: 0.2051

There is an increase of rainfree days between 2005 and 2022, but not statistically significant when zooming in to the shorter period of interest 2014-2021.

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 
-14.907  -7.157  -3.376   8.296  18.385 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) 2863.2250   957.7131   2.990  0.00866 **
years         -1.3540     0.4756  -2.847  0.01166 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.47 on 16 degrees of freedom
Multiple R-squared:  0.3362,    Adjusted R-squared:  0.2947 
F-statistic: 8.103 on 1 and 16 DF,  p-value: 0.01166

lm.weakrain.2014_2021 <- lm(weakraindays[yrs_2014_2018]~years[yrs_2014_2018])
summary(lm.weakrain.2014_2021)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-17.7143  -3.0982   0.9464   4.3036  11.6429 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)  
(Intercept)          6836.357   3135.137   2.181   0.0720 .
years[yrs_2014_2018]   -3.321      1.554  -2.137   0.0764 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.07 on 6 degrees of freedom
Multiple R-squared:  0.4323,    Adjusted R-squared:  0.3376 
F-statistic: 4.568 on 1 and 6 DF,  p-value: 0.07643

The number of days with weak rain decreases both over the long period and still a trend remains over the short period.

Substantial rain

plot(years, strongraindays)

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

Call:
lm(formula = strongraindays ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.9842 -4.7378 -0.8939  5.0718 12.0230 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  45.489852 636.611373   0.071    0.944
years        -0.007224   0.316170  -0.023    0.982

Residual standard error: 6.959 on 16 degrees of freedom
Multiple R-squared:  3.263e-05, Adjusted R-squared:  -0.06247 
F-statistic: 0.000522 on 1 and 16 DF,  p-value: 0.9821

The number of days with substantial precipitation did not change over the years. Note that this is not necessarily saying much, because strong rain can mean anything between 5 mm - 45 mm rainfall per day.

Note however, that the annual precipitation depends a great deal on how many days it rained above 5 mm/day and is also strongly determined by the total precpitation received during those days:

par(mfrow=c(1,2))
plot(P, P_strong/P ,
     xlab = "annual precipitation in mm/a",
     ylab = "total P on days P>5mm / total P")
plot(P,strongraindays/raindays,
     xlab = "annual precipitation in mm/a",
     ylab = "no of days with P>5mm / no precipitation days")

Only 12-24% of the precipitation days have precipitation > 5 mm, but they are responsible for 45 to 70% of the annual precipitation. If they are single events, they also roughly correspond to the event size required to yield any soil moisture response (see Fischer et al. 2023, in print in HESS). Their frequency does not change.

We can use the relation above to calculate the expected strong precipitation each year based on the annual precipitation. This allows a more focussed look on how the strong precipitation changes:

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

lm.Pstrong_P <- lm(P_strong~P) # Relation between total annual precipitation and strong events
summary(lm.Pstrong_P)

Call:
lm(formula = P_strong ~ P)

Residuals:
    Min      1Q  Median      3Q     Max 
-59.669 -16.413   3.817  19.657  30.012 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -136.37471   37.48263  -3.638  0.00221 ** 
P              0.89074    0.07368  12.089 1.85e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 24.56 on 16 degrees of freedom
Multiple R-squared:  0.9013,    Adjusted R-squared:  0.8951 
F-statistic: 146.1 on 1 and 16 DF,  p-value: 1.849e-09

# predict contribution of strong events from annual precipitation as a baseline.
P_strong_lm <- predict(lm.Pstrong_P) # this is the expected contribution of strong events to total P

d.pred.P_strog <- P_strong-P_strong_lm # prediction error of contribution of strong precip days

# check whether there is a trend in over/underprediction
plot(years, d.pred.P_strog,
     ylab = "observed-predicted")

lm.pred_trend <- lm(d.pred.P_strog ~ years)
summary(lm.pred_trend)

Call:
lm(formula = d.pred.P_strog ~ years)

Residuals:
    Min      1Q  Median      3Q     Max 
-47.959 -12.293   4.658  17.718  25.987 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)  
(Intercept) -4286.9821  1974.9825  -2.171   0.0454 *
years           2.1291     0.9809   2.171   0.0454 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 21.59 on 16 degrees of freedom
Multiple R-squared:  0.2275,    Adjusted R-squared:  0.1792 
F-statistic: 4.712 on 1 and 16 DF,  p-value: 0.04535

lm.pred_trend.2014_2021 <- lm(d.pred.P_strog[yrs_2014_2021] ~ years[yrs_2014_2021])
summary(lm.pred_trend.2014_2021)

Call:
lm(formula = d.pred.P_strog[yrs_2014_2021] ~ years[yrs_2014_2021])

Residuals:
    Min      1Q  Median      3Q     Max 
-35.600 -11.517   4.437  16.273  20.091 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)          -4593.657   6621.710  -0.694    0.514
years[yrs_2014_2021]     2.281      3.282   0.695    0.513

Residual standard error: 21.27 on 6 degrees of freedom
Multiple R-squared:  0.07452,   Adjusted R-squared:  -0.07973 
F-statistic: 0.4831 on 1 and 6 DF,  p-value: 0.513

There is a relation in the prediction error with time. In other words: The contribution of large rainfall events to total rainfall in recent years is higher than expected based on the entire sample (2005-2022). This relation is not at all obvious for the shorter period 2014-2021.