Loading [MathJax]/extensions/Safe.js

Vanessa Riegler and Simon Widy

Group project: the Climate System

Instructions

Objectives

In this final project, you will apply the methods you learned over the past weeks to answer the questions below.

Deadline

Please submit your project via OLAT before Thursday January 11 at 00H (in the night from Wednesday to Thursday).

Formal requirements

You will work in groups of two. If we are an odd number of students, one group can have three participants. (Tip: I recommend that students who have not followed a programming class to team up with students who have).

Each group will submit one (executed) jupyter notebook containing the code, plots, and answers to the questions (text in the markdown format). Please also submit an HTML version of the notebook. Each group member must contribute to the notebook. The notebook should be self-contained and the answers must be well structured. The plots must be as understandable as possible (title, units, x and y axis labels, appropriate colors and levels…).

Please be concise in your answer. We expect a few sentences per answer at most - there is no need to write a new text book in this project! Use links and references to the literature or your class slides where appropriate.

Grading

We will give one grade per project, according to the following table (total 10 points):

• correctness of the code and the plots: content, legends, colors, units, etc. (3 points)
• quality of the answers: correctness, preciseness, appropriate use of links and references to literature or external resources (5 points)
• originality and quality of the open research question (2 points)

# Import the tools we are going to need today:
import matplotlib.pyplot as plt  # plotting library
import numpy as np  # numerical library
import xarray as xr  # netCDF library
import pandas as pd  # tabular library
import cartopy  # Map projections libary
import cartopy.crs as ccrs  # Projections list
import cartopy.feature as cfeature
# Some defaults:
plt.rcParams['figure.figsize'] = (12, 5)  # Default plot size

Part 1 - temperature climatology

Open the ERA5 temperature data:

ds = xr.open_dataset('ERA5_LowRes_Monthly_t2m.nc')
ds['t2m'] -= 273.15 
ds['t2m'].attrs['units'] = 'degC'
ds = ds.rename({'t2m': 'mean_temp'})

Plot three global maps:

• Compute and plot the temporal mean temperature $\overline{T}$ for the entire period (unit °C)
• Compute and plot $\overline{T^{*}}$ (see lesson), the zonal anomaly map of average temperature.
• Compute and plot the monthly average temperature for each month $\overline{T_M}$ (annual cycle). I expect a variable of dimensions (month: 12, latitude: 241, longitude: 480). Hint: remember the .groupby() command we learned in the lesson. Now plot the average monthly temperature range, i.e. $\overline{T_M}max$ - $\overline{T_M}min$ on a map.

Questions:

1. Look at the zonal temperature anomaly map.
   • Explain why norther Europe and the North Atlantic region is significantly warmer than the same latitudes in North America or Russia.
   • Explain why the Northern Pacific Ocean does not have a similar pattern.
2. Look at the average monthly temperature range map.
   • Explain why the temperature range is smaller in the tropics than at higher latitudes
   • Explain why the temperature range is smaller over oceans than over land
   • Where is the temperature range largest? Explain why.

Answers:

1. Ocean currents:

• The presence of the Gulf Stream in the North Atlantic plays a crucial role. The Gulf Stream carries warm water from the tropics towards higher latitudes, particularly affecting the western coast of Europe. This warm ocean current helps moderate temperatures in Northern Europe.

• While the Gulf Stream transports warm waters into the North Atlantic, the North Pacific lacks an equivalent strong, warm ocean current reaching high latitudes. The Kuroshio Current in the North Pacific, although warm, tends to move along lower latitudes and has less influence on the temperature at higher northern latitudes compared to the Gulf Stream.

2. Average Monthly Temperature Range Map:

• Solar insolation: The tropics receive relatively consistent and high solar insolation throughout the year. The sun is nearly overhead, leading to consistent day length and more uniform heating.

• High Humidity: Tropical regions often have high humidity levels. Water vapor acts as a greenhouse gas, trapping heat and reducing temperature extremes.

• High specific heat of water: Water has a higher specific heat capacity than land. Oceans can absorb and store large amounts of heat without a significant increase in temperature. This moderating effect results in smaller temperature ranges over oceans.

• The temperature range is largest in [Sibirea](https://en.wikipedia.org/wiki/Siberia#:~:text=By%20far%20the%20most%20commonly,average%20about%2010%20%C2%B0C%20(), the combination of a continental climate, distance from oceans, low specific heat of land, Arctic influences, snow cover, and the region's latitude collectively result in the large temperature range observed in Siberia.

# Compute the temporal mean temperature
T_mean = ds['mean_temp'].mean(dim='time')

# Plot the temporal mean temperature
T_mean.plot()
plt.title('temporal mean temperature')
plt.show()

# Compute the zonal anomaly map of average temperature
T_anomaly = ds['mean_temp'].mean(dim='time') - ds['mean_temp'].mean(dim=['time', 'longitude'])

# Plot the zonal anomaly map of average temperature
T_anomaly.plot()
plt.title('zonal anomaly map of average temperature')
plt.show()

# Compute the monthly average temperature for each month
T_m = ds['mean_temp'].groupby('time.month').mean(dim='time')

# Plot annual cycle
fig, ax = plt.subplots(figsize=(12, 5))

T_m.mean(dim=['longitude', 'latitude']).plot(ax=ax)
ax.set_title('annual cycle')
ax.set_xlabel('month')
ax.set_ylabel('average temperature(°C)')

ax.set_xticks(np.arange(1, 13, 1))
ax.set_xticklabels(np.arange(1, 13, 1))
ax.set_yticks(np.arange(3, 9, 1))
ax.set_yticklabels(np.arange(3, 9, 1))

ax.grid(True)
plt.show()

# Compute the average monthly temperature range
T_m_range = T_m.max(dim='month') - T_m.min(dim='month')

# Plot the average monthly temperature range
T_m_range.plot(cmap='coolwarm')
plt.title('average monthly temperature range')
plt.show()

Part 2 - Precipitation climatology

Open the precipitation file and explore it. The units of monthly precipitation are wrongly labeled (unfortunately). They should read: m per day.

ds_prec = xr.open_dataset('ERA5_LowRes_Monthly_tp.nc')
ds_prec['tp'] *= 1000  # Convert to mm
ds_prec['tp'].attrs['units'] = 'mm'

Using .groupby(), compute the average daily precipitation for each month of the year (I expect a variable of dimensions (month: 12, latitude: 241, longitude: 480)). Convert the units to mm per day. Plot a map of average daily precipitation in January and in August with the levels [0.5, 1, 2, 3, 4, 5, 7, 10, 15, 20, 40] and the colormap `YlGnBu'

Questions:

1. Describe the location of the ITCZ in January and February. Without going into the details, explain (in one or two sentences)
2. Describe the precipitation seasonality in West Africa and in India. Name the phenomenon at play.

Answers:

1. In January and February, the [Intertropical Convergence Zone (ITCZ)](https://www.dwd.de/DE/service/lexikon/Functions/glossar.html?lv3=101278&lv2=101224#:~:text=F%C3%BCr%20ITCZ%20(englisch%3A%20Inter%2D,und%20kurze%20tropischen%20Gewitterst%C3%BCrme%20abwechseln.) generally lies near or slightly south of the equator, with the tilt of the Earth playing a role, especially in the Southern Hemisphere.

2. In West Africa, there's a distinct wet season (June to September) and dry season. The West African Monsoon is the driving phenomenon. In India, the Indian Monsoon brings the majority of rainfall from June to September, with a secondary contribution from the northeast monsoon (October to December).

month_dict = {
    1: 'January',
    2: 'February',
    3: 'March',
    4: 'April',
    5: 'May',
    6: 'June',
    7: 'July',
    8: 'August',
    9: 'September',
    10: 'October',
    11: 'November',
    12: 'December'
}


# Compute the average daily precipitation for each month
m_precip = ds_prec['tp'].groupby('time.month').mean('time')

levels = [0.5, 1, 2, 3, 4, 5, 7, 10, 15, 20, 40]
cmap = 'YlGnBu'

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12,10), subplot_kw={'projection':ccrs.PlateCarree()})

# Plot average daily precipitation for January
month = 1
m_precip.sel(month=month).plot(ax=ax1, levels=levels, cmap=cmap, transform=ccrs.PlateCarree())
ax1.coastlines()
gridlines = ax1.gridlines(draw_labels=True)
gridlines.bottom_labels  = False
gridlines.right_labels = False
ax1.set_title(f'average daily precipitation in {month_dict[month]}', rotation=90, x=-0.1, y=0)

# Plot average daily precipitation for August
month = 8
m_precip.sel(month=month).plot(ax=ax2, levels=levels, cmap=cmap, transform=ccrs.PlateCarree())
ax2.coastlines()
gridlines = ax2.gridlines(draw_labels=True)
gridlines.top_labels  = False
gridlines.right_labels = False
ax2.set_title(f'average daily precipitation in {month_dict[month]}', rotation=90, x=-0.1, y=0)

plt.show()

Part 3: sea-level pressure and surface winds

Open the file containing the surface winds (u10 and v10) and sea-level pressure (msl).

ds_wind = xr.open_dataset('ERA5_LowRes_Monthly_uvslp.nc')

Compute $\left[ \overline{SLP} \right]$ (the temporal and zonal average of sea-level pressure). Convert it to hPa, and plot it (line plot). With the help of plt.axhline, add the standard atmosphere pressure line to the plot to emphasize high and low pressure regions. Repeat with $\left[ \overline{u_{10}} \right]$ and $\left[ \overline{v_{10}} \right]$ (in m s$^{-1}$) and add the 0 horizontal line to the plot (to detect surface westerlies from easterlies for example).

Questions:

1. Based on your knowledge about the general circulation of the atmosphere, explain the latitude location of the climatological high and low pressure systems of Earth.
2. Similarly, explain the direction and strength of the zonal and meridional winds on earth (tip: the sea-level pressure plot helps)

Answers:

1. The general circulation of the Earth's atmosphere includes the equatorial low-pressure zone near the equator, subtropical high-pressure zones around 20-30 degrees latitude, and polar low-pressure zones near the poles. These patterns are attributed to the distribution of solar radiation.

2. Zonal winds (east-west) on Earth are influenced by pressure differences between the equator and mid-latitudes, creating Trade Winds near the equator and prevailing westerlies in mid-latitudes. Meridional winds (north-south) result from the Hadley and Ferrel Cell circulations, causing air to move poleward and equatorward. Near the poles, polar easterlies flow equatorward. These wind patterns are driven by the distribution of sea-level pressure, influenced by solar heating, and are affected by the Coriolis effect due to Earth's rotation.

# Compute the temporal and zonal average of sea-level pressure
a_slp = ds_wind['msl'].mean(dim=['time', 'longitude'])

# Convert it to hPa
a_slp_hpa = a_slp / 100

fig, ax1 = plt.subplots(figsize=(10, 6))

# Plot sea-level pressure
color = 'tab:red'
ax1.axhline(y=1013.25, color='r', linestyle='--', label='Standard Atmosphere Pressure')
ax1.set_xlabel('Latitude')
ax1.set_ylabel('Pressure (hPa)', color=color)
ax1.plot(a_slp_hpa['latitude'], a_slp_hpa, color=color, label='Sea Level Pressure')
ax1.tick_params(axis='y', labelcolor=color)

ax2 = ax1.twinx()

# Plot wind speeds
colors = ['blue', 'lightblue']
variables = ['u10', 'v10']
ax2.set_ylabel('Speed (m/s)')
ax2.axhline(y=0, color='b', linestyle='--', linewidth=1, label='0-line')

for var, color in zip(variables, colors):
    a_var = ds_wind[var].mean(dim=['time', 'longitude'])
    ax2.plot(a_var['latitude'], a_var, color=color, label=f'{var} Wind Speed')

ax1.set_ylim(bottom=a_slp_hpa.min(), top=a_slp_hpa.max())
ax2.set_ylim(bottom=-8, top=8)

ax1.grid(True)
ax2.grid(True)
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')

plt.title('temporal and zonal average of sea-level pressure and wind speeds')

fig.tight_layout()
plt.show()

Part 4: temperature change and CO$_2$ concentrations

Download the global average CO$_2$ concentration timeseries data in the CSV format (source: NOAA). Here, let me help your read them using pandas:

# Open csv as dataframe with matching length
df = pd.read_csv('co2_mm_gl.csv', skiprows=38, skipfooter=57, parse_dates={'date': [0, 1]}, index_col='date', engine='python', date_format='%Y %m')

# Add needed column to dataset
ds['co2'] = ('time', df['average'])

Prepare three plots:

• plot the monthly global CO$_2$ concentration as a function of time.
• plot the annual average timeseries of global CO$_2$ concentration as a function of time.
• plot the annual average timeseries of global CO$_2$ concentration and of global 2m temperature from ERA5 on the same plot (using a secondary y axis for temperature).

Questions:

1. Describe and explain the annual cycle of CO$_2$ concentrations
2. What was the CO$_2$ concentration in the atmosphere in the pre-industrial era? Compute the annual increase in CO$_2$ concentration (unit: ppm per year) between 1980 and 1985 and between 2016 and 2021.
3. Describe the relationship between global temperatures and CO$_2$ concentrations. Beside CO$_2$, name three processes that can influence temperature variability and change at the global scale.

Answers:

1. Annual cycle
2. Pre-industrial era

• Before the Industrial Revolution started in the mid-1700s, atmospheric carbon dioxide was 280 ppm or less.
• annual increase 1980-1985: 1.33ppm per year
• annual increase 2016-2021: 2.19ppm per year

3. • Relationship between global temperatures and CO2

• three processes:

• Solar Variability: Changes in solar radiation, such as variations in the sun's output, can influence global temperatures. While solar activity does play a role in climate variations, it is not the primary driver of recent warming trends observed on Earth.

• Volcanic Activity: Large volcanic eruptions release substantial amounts of aerosols and gases into the atmosphere. These particles can reflect sunlight back into space, leading to a temporary cooling effect known as volcanic cooling. However, the influence of volcanic activity on long-term temperature trends is relatively short-lived.

• Oceanic Processes: Oceanic circulation patterns, such as El Niño and La Niña events, can have a significant impact on global temperatures. These phenomena involve the periodic warming (El Niño) or cooling (La Niña) of sea surface temperatures in the central and eastern equatorial Pacific Ocean. They can influence weather patterns and temperature anomalies worldwide.

# Compute the annual increase in CO2 concentration (unit: ppm per year) between 1980 and 1985 and 2016-2021
df['annual_increase'] = df['average'].diff(periods=12)

filtered_df = df['1980-01-01':'1985-12-01']
filtered_df2 = df['2016-01-01':'2021-12-01']

total_increase = filtered_df['average'].iloc[-1] - filtered_df['average'].iloc[0]
total_increase2 = filtered_df2['average'].iloc[-1] - filtered_df2['average'].iloc[0]

average_annual_increase = total_increase / (filtered_df.index.year.max() - filtered_df.index.year.min() + 1)
average_annual_increase2 = total_increase2 / (filtered_df2.index.year.max() - filtered_df2.index.year.min() + 1)

print(f"The annual increase in CO2 concentration between 1980-1985: {average_annual_increase:.2f} ppm per year")
print(f"The annual increase in CO2 concentration between 2016-2021: {average_annual_increase2:.2f} ppm per year")

The annual increase in CO2 concentration between 1980-1985: 1.33 ppm per year
The annual increase in CO2 concentration between 2016-2021: 2.19 ppm per year

annual_a = ds['co2'].resample(time='Y').mean()

# Compute the global temp for every timestep resample by years and get means
a_temp = ds['mean_temp'].mean(dim=['longitude', 'latitude']).resample(time='Y').mean()

fig, ax1 = plt.subplots(figsize=(10, 6))

# Plot monthly and yearly average of CO2 concentration
color = 'tab:blue'
ax1.set_xlabel('Year')
ax1.set_ylabel('CO2 Concentration (ppm)', color=color)
ln1 = ax1.plot(ds['time'], ds['co2'], linewidth=1, color='black', label='Monthly CO2 Concentration')
ln2 = ax1.plot(annual_a['time'], annual_a, linewidth=3, color=color, label='Annual Mean CO2')
ax1.tick_params(axis='y', labelcolor=color)

# Create second y-axis
ax2 = ax1.twinx()

# Plot yearly average temperature and convert it from Kelvin to Celsius
color = 'tab:red'
ax2.set_ylabel('Temperature (°C)', color=color)
ln3 = ax2.plot(a_temp['time'], a_temp, linewidth=3, color=color, label='Annual Mean Temperature')
ax2.tick_params(axis='y', labelcolor=color)

plt.title('annual average global CO2 concentration and 2m temperature')

plt.grid(True)

lns = ln1 + ln2 + ln3
labs = [l.get_label() for l in lns]
ax1.legend(lns, labs, loc='upper left')

plt.show()

Part 5: variability and ENSO (open research question)

As of November 2023, the world is currently experiencing El Niño conditions, with projections indicating its persistence through the upcoming spring (source). To understand the implications and compare with La Niña, we will utilize our available data.

Using your learned skills, you should create global anomalie maps of sea-surface temperature, precipitation, and air temperature, comparing one El Niño with one La Niña year. Describe the anomaly patterns you are seeing in your plots (e.g. where do we see an increase/decrease in precipitation).

We suggest to look for literature or a google search on strong ENSO events in the past 40 years, and select one good example for a positive and a negative phase. With citation or links, explain why you picked these years as examples.

Answer:

Our example for a strong positive phase: El Niño 2009/2010 (https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2011GL048521)

Example for a strong negative phase: La Niña 2010-2012 (https://en.wikipedia.org/wiki/2010%E2%80%932012_La_Ni%C3%B1a_event) (We did not use data from 2010 for the La Niña event but exclusively focused on 2011-2012 to highlight the differences.)

Sea suface temperature:

La Niña 2010-2012: La Niña is characterized by cooler sea surface temperatures in the central and eastern Pacific. Cooler sea surface temperatures in the eastern Pacific lead to enhanced trade winds and increased upwelling of cool, nutrient-rich deep water.

The western parts of the Pacific, including the west coast of Australia, tend to be warmer as the trade winds push warm surface water westward.

El Niño 2009-2010: El Niño is characterized by warmer sea surface temperatures in the central and eastern Pacific. Increased sea surface temperatures result in altered atmospheric conditions, including weakened trade winds and a decrease in the upwelling of cool deep water.

During El Niño events, parts of the western Pacific, including the west coast of Australia, may experience cooler conditions.

In summary, La Niña tends to bring cooler sea surface temperatures to the central and eastern Pacific, while El Niño leads to warmer conditions.

Precipitation:

La Niña 2010-2012: During La Niña, there is typically above-average rainfall in the western Pacific, Southeast Asia, Australia, and parts of South America. Conversely, La Niña often leads to drier-than-normal conditions in the central and eastern equatorial Pacific region, as well as in some parts of North America.

El Niño 2009-2010: El Niño tends to bring above-average rainfall to the central and eastern equatorial Pacific region, parts of South America (including Peru and Ecuador), and the southern United States. Conversely, El Niño often leads to drier-than-normal conditions in Australia, Southeast Asia, and some regions of Africa.

Air temperature:

La Niña: La Niña tends to have a cooling effect on global temperatures. During La Niña, cooler-than-average sea surface temperatures in the central and eastern Pacific contribute to a suppression of global atmospheric temperatures. In North America, La Niña typically brings colder and drier conditions, but for example in parts of Northern Russia, may experience milder winters due to a strengthened low-pressure influence over the Arctic.

El Niño: El Niño is associated with a warming effect on global temperatures. During El Niño, warmer-than-average sea surface temperatures in the central and eastern Pacific contribute to an overall warming of the atmosphere. Regions like North America often experience unusually high temperatures, and some areas may suffer from drought conditions. Africa is often affected by severe drought and dryness. Concurrently, parts of Russia, especially in Siberia, might experience colder-than-average winters due to changes in atmospheric circulation patterns.

Summary: The 09/10 El Niño possessed a unique characteristic among other El Niño events in that it not only showed the record-breaking warm SST anomaly in the central Pacific but also showed the fast phase transition to La Niña. Due to that, we chose these two events and compared it, because it was really intresting to see, the fast changes in global atmospheric patterns.

#data: precipitation, air temperature and sea-surface temperature
ds_sst = xr.open_dataset('ERA5_LowRes_Monthly_sst.nc')
ds_sst['sst'] -= 273.15 
ds_sst['sst'].attrs['units'] = 'degC'

# El Niño 2009/10
el_nino_sst = ds_sst['sst'].sel(time=slice('2009-01-01', '2010-12-31'))
el_nino_prec = ds_prec['tp'].sel(time=slice('2009-01-01', '2010-12-31'))
el_nino_mean_temp = ds['mean_temp'].sel(time=slice('2009-01-01', '2010-12-31'))

# La Niña 2011/2012
la_nina_sst = ds_sst['sst'].sel(time=slice('2011-01-01', '2012-12-31'))
la_nina_prec = ds_prec['tp'].sel(time=slice('2011-01-01', '2012-12-31'))
la_nina_mean_temp = ds['mean_temp'].sel(time=slice('2011-01-01', '2012-12-31'))

# Anomalies
el_nino_sst_anomaly = el_nino_sst.mean(dim='time')-ds_sst['sst'].mean(dim='time')
la_nina_sst_anomaly = la_nina_sst.mean(dim='time')-ds_sst['sst'].mean(dim='time')

el_nino_prec_anomaly = el_nino_prec.mean(dim='time')-ds_prec['tp'].mean(dim='time')
la_nina_prec_anomaly = la_nina_prec.mean(dim='time')-ds_prec['tp'].mean(dim='time')

el_nino_mean_temp_anomaly = el_nino_mean_temp.mean(dim='time')-ds['mean_temp'].mean(dim='time')
la_nina_mean_temp_anomaly = la_nina_mean_temp.mean(dim='time')-ds['mean_temp'].mean(dim='time')

# Calculation of zonal means
zonal_mean_el_nino_sst = el_nino_sst.mean(dim='longitude')
zonal_mean_la_nina_sst = la_nina_sst.mean(dim='longitude')

zonal_mean_el_nino_prec = el_nino_prec.mean(dim='longitude')
zonal_mean_la_nina_prec = la_nina_prec.mean(dim='longitude')

zonal_mean_el_nino_mean_temp = el_nino_mean_temp.mean(dim='longitude')
zonal_mean_la_nina_mean_temp = la_nina_mean_temp.mean(dim='longitude')

# Calculation of zonal anomalies
zonal_anomaly_el_nino_sst = zonal_mean_el_nino_sst - zonal_mean_el_nino_sst.mean(dim='time')
zonal_anomaly_la_nina_sst = zonal_mean_la_nina_sst - zonal_mean_la_nina_sst.mean(dim='time')

zonal_anomaly_el_nino_prec = zonal_mean_el_nino_prec - zonal_mean_el_nino_prec.mean(dim='time')
zonal_anomaly_la_nina_prec = zonal_mean_la_nina_prec - zonal_mean_la_nina_prec.mean(dim='time')

zonal_anomaly_el_nino_mean_temp = zonal_mean_el_nino_mean_temp - zonal_mean_el_nino_mean_temp.mean(dim='time')
zonal_anomaly_la_nina_mean_temp = zonal_mean_la_nina_mean_temp - zonal_mean_la_nina_mean_temp.mean(dim='time')

# Plot
fig, axes = plt.subplots(3, 2, figsize=(22, 24), subplot_kw={'projection': ccrs.PlateCarree()})

# Define custom colormap
cmap = plt.cm.viridis
custom_cmap = cmap

# SST-Anomalies
contour_el_nino_sst = axes[0, 0].contourf(ds_sst['longitude'], ds_sst['latitude'], el_nino_sst_anomaly, transform=ccrs.PlateCarree(), cmap=custom_cmap, levels=15, extend='both')
axes[0, 0].set_title('SST Anomalies - El Niño 2009/10')
axes[0, 0].add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':')
axes[0, 0].set_xlabel('Longitude')
axes[0, 0].set_ylabel('Latitude')

contour_la_nina_sst = axes[0, 1].contourf(ds_sst['longitude'], ds_sst['latitude'], la_nina_sst_anomaly, transform=ccrs.PlateCarree(), cmap=custom_cmap, levels=15, extend='both')
axes[0, 1].set_title('SST Anomalies - La Niña 2011/2012')
axes[0, 1].add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':')
axes[0, 1].set_xlabel('Longitude')
axes[0, 1].set_ylabel('Latitude')

# Precipitation-Anomalies
contour_el_nino_prec = axes[1, 0].contourf(ds_prec['longitude'], ds_prec['latitude'], el_nino_prec_anomaly, transform=ccrs.PlateCarree(), cmap=custom_cmap, levels=15, extend='both')
axes[1, 0].set_title('Precipitation Anomalies - El Niño 2009/10')
axes[1, 0].add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':')
axes[1, 0].set_xlabel('Longitude')
axes[1, 0].set_ylabel('Latitude')

contour_la_nina_prec = axes[1, 1].contourf(ds_prec['longitude'], ds_prec['latitude'], la_nina_prec_anomaly, transform=ccrs.PlateCarree(), cmap=custom_cmap, levels=15, extend='both')
axes[1, 1].set_title('Precipitation Anomalies - La Niña 2011/2012')
axes[1, 1].add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':')
axes[1, 1].set_xlabel('Longitude')
axes[1, 1].set_ylabel('Latitude')

# Air Temperature-Anomalies
contour_el_nino_mean_temp = axes[2, 0].contourf(ds['longitude'], ds['latitude'], el_nino_mean_temp_anomaly, transform=ccrs.PlateCarree(), cmap=custom_cmap, levels=15, extend='both')
axes[2, 0].set_title('Air Temperature Anomalies - El Niño 2009/10')
axes[2, 0].add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':')
axes[2, 0].set_xlabel('Longitude')
axes[2, 0].set_ylabel('Latitude')

contour_la_nina_mean_temp = axes[2, 1].contourf(ds['longitude'], ds['latitude'], la_nina_mean_temp_anomaly, transform=ccrs.PlateCarree(), cmap=custom_cmap, levels=15, extend='both')
axes[2, 1].set_title('Air Temperature Anomalies - La Niña 2011/2012')
axes[2, 1].add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':')
axes[2, 1].set_xlabel('Longitude')
axes[2, 1].set_ylabel('Latitude')

# Add colorbars
cbar_sst = plt.colorbar(contour_el_nino_sst, ax=axes[0, :], orientation='horizontal', pad=0.05, aspect=40, shrink=0.8)
cbar_prec = plt.colorbar(contour_el_nino_prec, ax=axes[1, :], orientation='horizontal', pad=0.05, aspect=40, shrink=0.8)
cbar_mean_temp = plt.colorbar(contour_el_nino_mean_temp, ax=axes[2, :], orientation='horizontal', pad=0.05, aspect=40, shrink=0.8)


# Adjust layout for colorbars
plt.subplots_adjust(left=0.1, right=0.9, bottom=0.2, top=0.9, wspace=0.2, hspace=0.3)

plt.show()