[This lecture was partially created with the support of Artificial Intelligence tools (ChatGPT)]
Introduction
The Po River Basin, located in northern Italy, is the country’s largest and most important drainage area. It covers roughly 70,000 square kilometers and extends from the western Alps to the Adriatic Sea. The Po River itself is about 652 kilometers long, making it Italy’s longest river. Major tributaries include the Dora Baltea, Ticino, Adda, Oglio, and Mincio rivers. The basin encompasses several major cities, such as Turin, Milan, and Ferrara.
It is a highly fertile region, forming the core of Italy’s agricultural production—especially rice, maize, and wheat. Intensive irrigation and industrial development have made it an economic powerhouse. However, these activities have also caused pollution and significant water management challenges. The area is prone to both flooding and drought due to climate variability and human pressures. Today, sustainable management of the Po River Basin is crucial for Italy’s environmental and economic future.
Climate change is significantly altering rainfall patterns and water availability across the basin. Rising temperatures increase evaporation and reduce snow accumulation in the Alps, which traditionally feeds the river during warmer months. As a result, summer river flows are decreasing, leading to more frequent and severe droughts. One example is the extreme low water levels recorded during the 2022 European drought. Reduced water flow threatens irrigation systems that support intensive farming in the Po Valley. Climate change also increases the risk of floods due to more intense rainfall events. Saltwater intrusion from the Adriatic Sea is becoming more common when river discharge is low. These environmental pressures affect biodiversity, agriculture, and hydroelectric power generation in the basin. Adapting water management strategies is therefore essential to protect the long-term sustainability of the Po River Basin.
In particular, floods over the Po River Basin are widely studied. The Po River is intensively monitored; in particular, river flow is monitored in several river cross sections. At Pontelagoscuro, which is considered the closure of the catchment, daily river levels are monitored since the 19th century. Since 1920 mean daily levels are converted into mean daily flows. the resulting time series is available here.
Note that the mean daily flow is not equivalent to the peak flow that occurred that day, but is lower. For a large basin like the Po an empirical relationship used in Italy to give the ratio between mean daily flow and peak flow in the day of a large discharge is given by the relationship
\(\frac{Q_{m}}{Q_{c}} \approx 0.75 - 0.85\) ,
that is called "Formula di Cotecchia". For the Po we usually adopt a value around 0.75.
Estimation of the peak value of the mean daily flow for the Po River at Pontelagoscuro
The peak river flow is usually estimated by fitting a proper extreme value probability distribution. Let us estimate it by fitting a Gumbel, a Frechet and a GEV distribution through the method of annual maxima and evaluate the goodness of the fit. Let us start with R.
# Convert Date column to Date format
poflow$Date=as.Date(poflow$Date, format = "%d/%m/%Y")
#Extract annual maxima
poflow$Year=format(poflow$Date, "%Y")
annual_max=aggregate(Flow ~ Year, data = poflow, max)
Now let's install the required packages to fit the distributions.
install.packages(c("evd","fitdistrplus","ismev","goftest"))
library(evd)
library(fitdistrplus)
library(ismev)
library(goftest)
Let's estimate the Gumbel and GEV distributions.
gumbel_fit <- fgev(annual_max$Flow, shape = 0)
summary(gumbel_fit)
#Fit a GEV distribution
gev_fit <- fgev(annual_max$Flow)
summary(gev_fit)
For the Frechet we need to apply maximum likelihood with an optimisation routine to constrain shape parameter. Note: given that GEV estimation ended up with a negative estimate for shape, it is likely that Frechet converges to shape zero.
#We need to ensure that shape parameter is greater than zero. Then, we need a custom maximum likelihood estimation with a lower bound on the shape parameter
#We first fit a GEV because we need an output with a specific format to make diagostic with the loaded packages
frechet_fit <- fgev(annual_max$Flow)
#By the way, if the shape parameter is greater than zero then we can stop here, the estimated GEV is already a Frechet
#Define parameter vector
frechet_fit1=list()
frechet_fit1$param=c(mean(annual_max$Flow),sd(annual_max$Flow),0.1)
# Define the function giving the likelihood
negloglik <- function(par, x)
{
loc=par[1]
scale=par[2]
shape=par[3]
-sum(dgev(x, loc=loc, scale=scale, shape=shape, log=TRUE))
}
# Impose Lower bounds: shape > 0, scale > 0
frechet_fit1$param=optim(par=frechet_fit$param,fn=negloglik,x=annual_max$Flow,method="L-BFGS-B",lower=c(-Inf, 1e-6, 1e-6))$par
#Copy parameter values into the dummy results obtained with GEV
frechet_fit[1]$estimate=frechet_fit1$param
Let's check goodness of the fit with the Kolmogorov-Smirnov test and probability plots.
#Gumbel
ks.test(annual_max$Flow,"pgumbel",loc=gumbel_fit$estimate["loc"],scale = gumbel_fit$estimate["scale"])
#If you get a warning about ties, it is due to the presence of duplicate values in the series of annual maxima
#Frechet
ks.test(annual_max$Flow,"pfrechet",loc=frechet_fit$estimate[1],scale=frechet_fit$estimate[2],shape = frechet_fit$estimate[3])
#GEV
ks.test(annual_max$Flow,"pgev",loc=gev_fit$estimate["loc"],scale = gev_fit$estimate["scale"],shape = gev_fit$estimate["shape"])
#Diagnostic plots
par(mfrow=c(2,2))
plot(gumbel_fit)
plot(gev_fit)
plot(frechet_fit)
For the sake of comparison, let's compute the peak flow for a given return period
T <- c(30,50,100,200)
#Design floods (GEV)
Qgev <- qgev(1-1/T,loc=gev_fit$estimate["loc"],scale=gev_fit$estimate["scale"],shape=gev_fit$estimate["shape"])
#Design floods (Gumbel)
Qgumbel <- qgumbel(1-1/T,loc=gumbel_fit$estimate["loc"],scale=gumbel_fit$estimate["scale"])
# Results table
results <- data.frame(ReturnPeriod=T,Q_GEV=Qgev,Q_Gumbel=Qgumbel)
print(results)
To conclude, fit Gumbel with the method of moments and compare the resulting parameter values and goodness of fit tests. You may also compare with the "homemade" estimation of distribution parameters that you can find here.
Last modified on March 16, 2026
- 15 viste




