Nonstationary FFA

This vignette demonstrates how to use the ffaframework package to perform flood frequency analysis (FFA) under the assumption of nonstationarity. Readers unfamiliar with stationary FFA workflows should first consult the Stationary FFA vignette.

The framework supports three forms of nonstationarity, modelled as linear trends in the distribution parameters:

A linear trend in the mean (location parameter).
A linear trend in the variance (scale parameter).
A linear trend in both the mean and variance.

The shape parameter is treated as a constant considering its difficult estimation due to sample size limitations.

Case Study

This vignette will explore data from the Bow River at Banff (05BB001) hydrological monitoring station. The remoteness of this station means that trends annual maxima are caused by changes in climate as opposed to changes in land use or cover. Data for this station is provided as Application_3.1.csv in the ffaframework package.

library(ffaframework)

csv_path <- system.file("extdata", "Application_3.1.csv", package = "ffaframework")
df <- read.csv(csv_path, comment.char = "#")
df <- subset(df, !is.na(max)) # Remove missing values

head(df)
#>    year max
#> 10 1909 314
#> 11 1910 230
#> 12 1911 264
#> 13 1912 174
#> 14 1913 232
#> 15 1914 214

plot_ams_data(df$max, df$year, title = "Bow River at Banff (05BB001)")

The `trend` List

This vignette assumes a scenario where a trend in the mean has been identified. To specify a trend in the time series during FFA, create a trend list as shown below:

trend <- list(location = TRUE, scale = FALSE)

Note: For guidance on trend detection, refer to the Trend Identification vignette.

Distribution Selection

L-moment-based distribution selection remains applicable under nonstationarity, but requires dataset decomposition (detrending) prior to analysis. This is accomplished using the ams_decomposition function, which takes an annual maximum series (AMS), the corresponding vector of years, and the trend object. The decomposed vector of data is then passed to the selection function.

data_decomposed <- ams_decomposition(df$max, df$year, trend)

selection <- select_ldistance(data_decomposed)

print(selection$recommendation)
#> [1] "GNO"

plot_lmom_diagram(selection)

Conclusion: The generalized normal (GNO) distribution is the best-fit distribution for the sample.

Parameter Estimation

Because L-moments parameter estimation requires stationarity, we use maximum likelihood estimation for nonstationary models. The fit_maximum_likelihood function implements maximum likelihood estimation for both stationary and nonstationary distributions. It has two required arguments:

data: The annual maximum series observations.
model: A three-letter code corresponding to a probability distribution (ex. "GNO").

Since a nonstationary model is used, two additional arguments are required:

years: The corresponding vector of years for the observations in data.
trend: The nonstationary trend object described above.

fit <- fit_maximum_likelihood(
    df$max,
    "GNO",
    years = df$year,
    trend = trend
)

print(fit$params)
#> [1] 224.0496619 -35.5153685  54.6324886  -0.3689085

print(fit$mll)
#> [1] -590.7329

Note: The fitted parameters are: \((\mu_0, \mu_1, \sigma, \kappa)\), where the time-dependent location is modeled as \(\mu(t) = \mu_0 + \mu_1 t\). The covariate \(t\) is a linear function of the years: \((\text{years} - 1900) / 100\).

Uncertainty Quantification

Uncertainty quantification is also essential for nonstationary probability distributions. In addition to the parametric bootstrap, the framework implements the regula-falsi profile likelihood (RFPL) method for MLE. The uncertainty_rfpl method has two required arguments:

data: The annual maximum series observations.
model: A three-letter code corresponding to a probability distribution (ex. "GNO").

Since the dataset has a nonstationary trend, three additional arguments are required:

years: The corresponding vector of years for the observations in data.
trend: The nonstationary trend object described above.
slices: The years at which return levels are computed.

uncertainty <- uncertainty_rfpl(
    df$max,
    "GNO",
    years = df$year,
    trend = trend,
    slices = c(1925, 2025)
)

print(uncertainty[[2]]$estimates)
#> [1] 179.6555 233.5731 269.1684 303.2462 347.4806 380.9029

plot_nsffa(uncertainty)

Example Conclusion: In the year 2025, there is a \(1/20\) exceedance probability of a flood of approximately \(300\text{m}^3/\text{s}\) or greater.