Introduction

The recovery of the hole is progressing. The annual data on the area of the ozone hole (from 07 September–13 October of each year) and the minimum of the mean Southern Hemisphere ozone concentration (from 21 September–16 October of each year) since 1979 is available on the NASA Ozone Watch website, which is updated each year.

The ozone concentration data shows a fair bit of year-to-year variability, especially in recent years. It can therefore be useful to apply a line of best fit to this data to get a better idea of the trend over the period since the introduction of the Montreal Protocol.

The data

The data from comes from the NASA Ozone Watch site: ozone.csv. It contains the area of ozone hole in km2, and the minimum ozone concentration in DU, for the period between 1979 to present.

library(skimr)

ozone <- read.csv("data/ozone.csv")

# skim(ozone)

You can see that …

Plotting the data

Minimum ozone concentration over time

plot(minimum ~ Year,
     data = ozone,
     type = "l")

This plot shows the clear decrease in ozone minimum from about 225 DU in 1979 to around 100 DU in the mid to late 1990s, followed by a gradual recovery. However, even up to 2020 there are years where the minimum ozone is still below 120 DU, and the values since about 2000 are quite variable, spanning a range of about 60 DU. This year-to-year variability can make it difficult to see the underlying trend clearly.

Ozone hole area over time

plot(area ~ Year,
     data = ozone,
     type = "l")

Rate of change since the introduction of the Montreal Protocol

It is evident that the introduction of the Montreal Protocol has halted the spring degradation of the ozone layer in the Southern Hemisphere, but to what extent is it recovering so far? We can quantify the recovery of the ozone layer by fitting a line of best fit to the data since the introduction of the Protocol. The slope of this line tells us the rate of change in ozone hole parameters.

The Protocol was introduced in 1990, and CFCs were completely phased out in developed countries by 1996. We can see from graphs above the the trajectory of ozone started to change in the mid-1990s. So we will use 1996 as the starting point.

Recall: selecting a subset of rows from a data frame, e.g….

subset <- dataset[dataset$var <= 10, ]

We extract data for the years 1996 on-wards into another dataset called recovery:

recovery <- ozone[ozone$Year >= 1996,]

The next step is to plot the area and minimum data (separately) and fit a line of best fit to the data.

area.lm <- lm(area ~ Year, data = recovery)

plot(area ~ Year,
     data = recovery,
     type = 'l',
     xlab ='Year',
     ylab ='Ozone hole area (million km^2)')

abline(area.lm,
       col = "darkred")

We will fit a line for the relationship between area and Year, using the lm function, and add it to the plot using abline().

area.lm <- lm(area ~ Year, data = recovery)

The fitted line shows a gradual decrease over time. The slope of this line represents the average rate that the size of the ozone hole has changed (decreased) since 1996, in million km2 per year. This rate (slope) is given by the coefficient \(a\) in the equation for a line \(y = a.x + b\). We can extract the coefficients from an lm() object using $coefficients.

area.lm$coefficients
## (Intercept)        Year 
## 254.8248721  -0.1159288

The coefficient of Year is the slope of the line, which should be a negative value as the values of y (the area) are decreasing as x (year) is increasing. Although this value appears small recall that the units for the area are million km2. Therefore, at the time of writing (2025), the area of the ozone hole is decreasing by about 116,000 km2 per year. For reference, the area of Bulgaria is about 110,000 km2.

Reapitting this analysis for the minimum ozone concentration:

# line of best fit to minimum ozone concentration
minimum.lm <- lm(minimum ~ Year,data = recovery)

# plot minimum ozone concentration against year and add the line of best fit
plot(minimum ~ Year,
     data = recovery,
     type  = "l")

abline(minimum.lm, col = "darkred")

Now extract the rate that the minimum ozone concentration is changing from the fit. What is the rate of change in DU per decade?

minimum.lm$coefficients
## (Intercept)        Year 
## -1186.55796     0.64901

Summary

The above analysis shows that, over the 30 years or so since the introduction of the Montreal Protocol, the overall trend is a decrease in the size of ozone hole and and an increase in the annual minimum ozone concentration.

The aim is return to 1980 levels which, for the minimum ozone concentration, was 203 DU. The expectation is that this will be reached (for average October total column ozone) by about 2060-2070. Looking at the trend so far, how optimistic are you that this will be achieved? What are the factors and uncertainties affecting the rate of recovery?

Adding a smoothed line to the data

Sometimes it is helpful to add a ‘smoothed’ line to plots where the data is quite variable and where the trend is changing over time, such as the full time series of the ozone data above. An example of this is shown in Figure 2.4.20 with the ‘Mean of observations’ line. Here you will learn how to apply a smoothing function to the data. Applying these smoothing trends are useful for an initial exploration and visualisation of the longer-term pattern. One appraoch that can be done in R quite easily the LOESS (LOcal regreESSion) smoothing method. LOESS works by fitting simple polynomial regressions to small, overlapping segments of the data and then combining these local fits into a smooth curve.

It is important to note that LOESS smoothing does not test for the statistical significance of a trend. Instead, it provides a helpful first step for highlighting the underlying trajectory in noisy data, which is typical of environmental time series data, and for guiding further, more formal trend analysis.

In the following code chunk the minimum ozone time series is plotted again and a LOESS smoothing curve is applied and added to the plot. First, an object called smooth.min is created, which contains a LOESS model that fits the minimum ozone values as a function of Year. - next, the smoothed values are calculated using the predict() function applied to smooth.min and the smoothed estimates are added as a line to the plot inside the lines function.

plot(minimum ~ Year,
     data = ozone,
     type = 'l',
     lty = 2,
     xlab = 'Year',
     ylab = 'Minimum ozone (DU)')

smooth.min <- loess(minimum ~ Year, data = ozone)

lines(ozone$Year, predict(smooth.min), lwd = 2, col = 'blue')

It is not essential that you understand the inner workings of these new functions. What matters is recognising that they create a smooth curve that passes through the central tendency of the data, helping us to visualise the overall trend despite the year-to-year variability.

Below is the same analysis for the area of the ozone hole. Call your LOESS model of the area data smooth.area.

plot(area ~ Year, data = recovery
     )

How do these smoothed trends compare to the linear trends calculate above?

plot(area ~ Year, data = ozone, type = 'l', lty=2, xlab='Year', ylab='Ozone hole area (million km^2)')
smooth.area <- loess(area ~ Year, data = ozone)
lines(ozone$Year, predict(smooth.area), lwd=2, col='blue')