Menghitung Lapse Rate dan Stabilitas dari ERA5

Sumber: NASA Earth Observatory / Jesse Allen / MODIS Rapid Response Team (halaman sumber)

Mengapa Stabilitas Atmosfer Penting

Indonesia berada di jantung Maritime Continent, salah satu kawasan dengan aktivitas konvektif tertinggi di dunia. Setiap sore, cumulonimbus tumbuh di atas Jawa, Sumatera, dan Kalimantan — menghasilkan hujan lebat, petir, dan banjir bandang yang kerap terjadi hanya dalam hitungan jam. Pertanyaan operasional yang dihadapi forecaster BMKG adalah: kapan atmosfer benar-benar siap "meledak" menjadi badai?

Jawabannya ada pada stabilitas atmosfer. Jika kolom udara berada dalam kondisi conditionally unstable, sebuah trigger kecil — orografi, konvergensi angin, atau pemanasan permukaan — sudah cukup untuk memulai konveksi dalam. Penelitian BMKG tentang Bandung Basin menunjukkan bahwa indeks-indeks stabilitas seperti CAPE dan Lifted Index menunjukkan tren naik 2–5 jam sebelum onset hujan lebat dan petir, menegaskan nilai diagnostik mereka untuk peringatan dini jangka pendek.

Tutorial ini menghitung environmental lapse rate (\(\Gamma\)) dari ERA5 cache, mengklasifikasikan setiap kolom ke tiga kelas stabilitas, dan membangun proxy buoyancy energy sederhana — tanpa sounding penuh atau API tambahan.

Lapse Rate, Adiabatic Lifting, dan Definisi Stabilitas

Environmental lapse rate (\(\Gamma\)) adalah laju perubahan suhu aktual atmosfer terhadap ketinggian, didefinisikan sebagai:

$$\Gamma = -\frac{dT}{dz}$$

Tanda negatif ada karena suhu umumnya turun seiring naiknya ketinggian, sehingga \(\Gamma\) bernilai positif. Kita mengukur \(\Gamma\) dari dua lapisan tekanan — 500 hPa dan 850 hPa — menggunakan data reanalysis ERA5.

Dua nilai referensi penting:

• Dry adiabatic lapse rate (DALR): \(\Gamma_d = g/c_p = 9{,}80665\,/\,1005 \approx 9{,}76\ \text{K/km}\) (umum dibulatkan ke \(9{,}8\)). Ini adalah laju pendinginan udara tak jenuh yang terangkat adiabatik — konstanta fisika.
• Moist adiabatic lapse rate (MALR): \(\Gamma_m \approx 6\ \text{K/km}\) sebagai aproksimasi mid-troposphere tropis. Nilainya lebih kecil dari DALR karena panas laten kondensasi sebagian mengimbangi pendinginan adiabatik. Perlu dicatat bahwa MALR sebenarnya bervariasi dengan suhu dan tekanan — sekitar 4 K/km di udara hangat tropis dan 9 K/km di udara dingin polar — sehingga nilai 6 K/km adalah aproksimasi tengah-tropis yang praktis.

Dari perbandingan \(\Gamma\) terhadap \(\Gamma_d\) dan \(\Gamma_m\), ada tiga kelas stabilitas:

Bagan klasifikasi stabilitas atmosfer berdasarkan perbandingan environmental lapse rate dengan \(\Gamma_d\) dan \(\Gamma_m\).

Wilayah Indonesia didominasi conditionally unstable karena udara lembap dari Samudra Hindia dan Samudra Pasifik selalu tersedia sebagai "bahan bakar" konveksi.

Membuka Data ERA5 dari Cache Lokal

Data ERA5 yang kita gunakan sudah tersimpan di cache lokal — tidak perlu koneksi ke CDS atau registrasi akun. File yang relevan adalah:

• /data/era5/era5_t_pl500-850_indonesia_2024_d.nc — suhu (K) pada 500 dan 850 hPa, harian 00 UTC
• /data/era5/era5_z_pl500-850_indonesia_2024_d.nc — geopotential (\(\text{m}^2/\text{s}^2\)) pada level yang sama

Keduanya mencakup bbox Indonesia [6°N, 95°E, −11°S, 141°E] untuk seluruh tahun 2024. Snippet di bawah membuka kedua file dan menampilkan struktur dataset sehingga kita bisa konfirmasi dimensi dan koordinat yang tersedia.

import xarray as xr
import numpy as np

# Buka kedua file ERA5 dari cache lokal
ds_t = xr.open_dataset("/data/era5/era5_t_pl500-850_indonesia_2024_d.nc")
ds_z = xr.open_dataset("/data/era5/era5_z_pl500-850_indonesia_2024_d.nc")

print("=== Dataset Suhu (T) ===")
print(f"Variabel  : {list(ds_t.data_vars)}")
print(f"Dimensi   : {dict(ds_t.dims)}")
print(f"Level     : {ds_t.pressure_level.values} hPa")
print(f"Rentang waktu: {str(ds_t.valid_time.values[0])[:10]} s/d {str(ds_t.valid_time.values[-1])[:10]} (UTC)")
print(f"Lat range : {float(ds_t.latitude.min()):.1f} — {float(ds_t.latitude.max()):.1f} °")
print(f"Lon range : {float(ds_t.longitude.min()):.1f} — {float(ds_t.longitude.max()):.1f} °")

print("\n=== Dataset Geopotential (Z) ===")
print(f"Variabel  : {list(ds_z.data_vars)}")
print(f"Dimensi   : {dict(ds_z.dims)}")
print(f"Level     : {ds_z.pressure_level.values} hPa")

=== Dataset Suhu (T) ===
Variabel  : ['t']
Dimensi   : {'valid_time': 366, 'pressure_level': 2, 'latitude': 69, 'longitude': 185}
Level     : [850. 500.] hPa
Rentang waktu: 2024-01-01 s/d 2024-12-31 (UTC)
Lat range : -11.0 — 6.0 °
Lon range : 95.0 — 141.0 °

=== Dataset Geopotential (Z) ===
Variabel  : ['z']
Dimensi   : {'valid_time': 366, 'pressure_level': 2, 'latitude': 69, 'longitude': 185}
Level     : [850. 500.] hPa

Output mengonfirmasi dua pressure level (500/850 hPa), 366 time step harian 00 UTC, dan grid bbox Indonesia, dengan variabel t (suhu) dan z (geopotential).

Menghitung Lapse Rate dari ERA5

Kita pilih tanggal 15 Januari 2024 00 UTC sebagai contoh — bulan Januari merupakan puncak musim hujan di sebagian besar Jawa dan Sumatera, sehingga diharapkan kondisi instabilitas dominan.

Langkah konversi unit yang perlu diperhatikan:

1. Suhu: ERA5 menyimpan \(T\) dalam Kelvin. Untuk lapse rate dalam K/km, kita bisa langsung pakai Kelvin karena selisih suhu tidak bergantung offset satuan (\(\Delta T_K = \Delta T_{°C}\)).
2. Geopotential ke geopotential height: \(h = \Phi / g\) dengan \(g = 9{,}80665\ \text{m/s}^2\). Hasilnya dalam meter (geopotential metres, gpm), yang secara praktis setara dengan meter geometrik di ketinggian ini.
3. Layer thickness: \(\Delta z = z_{500} - z_{850}\) dalam meter. Karena 500 hPa selalu lebih tinggi dari 850 hPa, \(\Delta z > 0\).
4. Lapse rate: \(\Gamma = -(T_{500} - T_{850}) / \Delta z \times 1000\) dalam K/km. Tanda negatif memastikan \(\Gamma\) positif ketika suhu menurun ke atas (kondisi normal).

# Konstanta
g = 9.80665  # m/s²
cp = 1005.0  # J/(kg·K)

# Pilih tanggal analisis
TARGET_DATE = "2024-01-15"

# Ekstrak T pada 500 dan 850 hPa (K)
t500 = ds_t["t"].sel(pressure_level=500, valid_time=TARGET_DATE, method="nearest")
t850 = ds_t["t"].sel(pressure_level=850, valid_time=TARGET_DATE, method="nearest")

# Ekstrak geopotential (m²/s²) dan konversi ke geopotential height (m)
z500_phi = ds_z["z"].sel(pressure_level=500, valid_time=TARGET_DATE, method="nearest")
z850_phi = ds_z["z"].sel(pressure_level=850, valid_time=TARGET_DATE, method="nearest")

z500_m = z500_phi / g  # geopotential height dalam meter
z850_m = z850_phi / g

# Layer thickness (m) — positif karena z500 > z850
delta_z = z500_m - z850_m

# Environmental lapse rate (K/km)
# Γ = -(T500 - T850) / (z500 - z850) * 1000
gamma = -(t500 - t850) / delta_z * 1000.0

# Statistik ringkas di atas grid Indonesia
gamma_vals = gamma.values.flatten()
gamma_vals = gamma_vals[~np.isnan(gamma_vals)]

print(f"Tanggal analisis    : {TARGET_DATE} 00 UTC")
print(f"Rata-rata Δz (500–850 hPa) : {float(delta_z.mean()):.0f} m")
print(f"\nDistribusi Environmental Lapse Rate Γ (K/km):")
print(f"  Min   : {gamma_vals.min():.2f}")
print(f"  P5    : {np.percentile(gamma_vals, 5):.2f}")
print(f"  P25   : {np.percentile(gamma_vals, 25):.2f}")
print(f"  Mean  : {gamma_vals.mean():.2f}")
print(f"  Median: {np.median(gamma_vals):.2f}")
print(f"  P75   : {np.percentile(gamma_vals, 75):.2f}")
print(f"  P95   : {np.percentile(gamma_vals, 95):.2f}")
print(f"  Max   : {gamma_vals.max():.2f}")
print(f"  Std   : {gamma_vals.std():.2f}")
print(f"\nCatatan: Γd = {g/cp*1000:.1f} K/km, Γm ≈ 6.0 K/km")

Tanggal analisis    : 2024-01-15 00 UTC
Rata-rata Δz (500–850 hPa) : 4375 m

Distribusi Environmental Lapse Rate Γ (K/km):
  Min   : 4.42
  P5    : 4.81
  P25   : 4.99
  Mean  : 5.11
  Median: 5.11
  P75   : 5.23
  P95   : 5.40
  Max   : 6.01
  Std   : 0.18

Catatan: Γd = 9.8 K/km, Γm ≈ 6.0 K/km

Dari statistik di atas, mean \(\Gamma \approx 5{,}11\ \text{K/km}\) — bahkan P95 hanya \(5{,}40\ \text{K/km}\) — berada di bawah threshold \(\Gamma_m = 6{,}0\ \text{K/km}\), sehingga praktis seluruh grid (99,99%) jatuh ke kelas absolutely stable. Ini bukan kesalahan data: snapshot harian 00 UTC menghaluskan boundary layer superadiabatik siang hari yang sesungguhnya memicu konveksi. Nilai \(\Delta z\) rata-rata \(4375\ \text{m}\) menjadi penyebut dalam kalkulasi \(\Gamma\) dan CAPE proxy, sehingga akurasi konversi geopotential penting.

Klasifikasi Stabilitas dan Integrasi CAPE Proxy

Dengan nilai \(\Gamma\) yang sudah terhitung di setiap grid point, kita lanjut ke dua tugas berikutnya: (1) mengklasifikasikan setiap kolom ke dalam tiga kelas stabilitas, dan (2) menghitung proxy buoyancy energy untuk mengukur potensi konvektif.

Untuk CAPE proxy, kita gunakan pendekatan dua-level yang sederhana. Sebuah parsel diangkat secara dry-adiabatik dari 850 hPa ke 500 hPa. Suhu parsel di 500 hPa adalah:

$$T_{\text{parcel,500}} = T_{850} - \Gamma_d \cdot \Delta z$$

Jika \(T_{\text{parcel,500}} > T_{500}\) (suhu parsel lebih hangat dari environment), parsel mengalami buoyancy positif. Energi proxy per kolom dihitung sebagai:

$$\text{CAPE proxy} = \frac{g}{T_{\text{avg}}} \cdot (T_{\text{parcel,500}} - T_{500}) \cdot \Delta z$$

dengan \(T_{\text{avg}} = (T_{500} + T_{850}) / 2\) sebagai suhu referensi lapisan, hasilnya dalam J/kg. CAPE operasional memerlukan sounding penuh dengan integrasi dari Level of Free Convection (LFC) hingga Equilibrium Level (EL); proxy dua-level ini hanya mengilustrasikan konsep buoyancy energy.

# Threshold klasifikasi (K/km)
GAMMA_M = 6.0   # moist adiabatic lapse rate (aproksimasi)
GAMMA_D = g / cp * 1000.0  # dry adiabatic: 9.80665/1005 * 1000 ≈ 9.76 K/km

# ---- Klasifikasi Stabilitas ----
stable = gamma < GAMMA_M
conditional = (gamma >= GAMMA_M) & (gamma <= GAMMA_D)
unstable = gamma > GAMMA_D

n_total = int(stable.count())
n_stable = int(stable.sum())
n_cond = int(conditional.sum())
n_unstab = int(unstable.sum())

print("=== Klasifikasi Stabilitas (2024-01-15 00 UTC) ===")
print(f"Total grid points       : {n_total}")
print(f"Absolutely stable       : {n_stable:5d}  ({100*n_stable/n_total:.1f}%)")
print(f"Conditionally unstable  : {n_cond:5d}  ({100*n_cond/n_total:.1f}%)")
print(f"Absolutely unstable     : {n_unstab:5d}  ({100*n_unstab/n_total:.1f}%)")
print(f"\nΓd digunakan = {GAMMA_D:.2f} K/km  |  Γm digunakan = {GAMMA_M:.1f} K/km")

# ---- CAPE Proxy (dua-level, dry adiabatic) ----
# Parsel diangkat dari 850 hPa ke 500 hPa secara dry-adiabatik
t_parcel_500 = t850 - (GAMMA_D / 1000.0) * delta_z  # K, lakukan dalam K/m lalu kalikan dz (m)

# Buoyancy temperature difference (K)
delta_T = t_parcel_500 - t500  # positif = buoyant

# Suhu rata-rata lapisan (K) sebagai referensi
t_avg = (t500 + t850) / 2.0

# CAPE proxy (J/kg) — integrasi dua-level
cape_proxy = (g / t_avg) * delta_T * delta_z

cp_vals = cape_proxy.values.flatten()
cp_vals = cp_vals[~np.isnan(cp_vals)]

print("\n=== CAPE Proxy Dua-Level (J/kg) ===")
print("  (Parsel dry-adiabatik dari 850 hPa, dibandingkan environment di 500 hPa)")
print(f"  Min   : {cp_vals.min():.1f}")
print(f"  P5    : {np.percentile(cp_vals, 5):.1f}")
print(f"  P25   : {np.percentile(cp_vals, 25):.1f}")
print(f"  Mean  : {cp_vals.mean():.1f}")
print(f"  Median: {np.median(cp_vals):.1f}")
print(f"  P75   : {np.percentile(cp_vals, 75):.1f}")
print(f"  P95   : {np.percentile(cp_vals, 95):.1f}")
print(f"  Max   : {cp_vals.max():.1f}")
print(f"\n  Fraksi kolom dengan CAPE proxy > 0 J/kg : {100*(cp_vals>0).mean():.1f}%")
print("  Catatan: CAPE operasional butuh sounding penuh + LFC/EL, bukan hanya 2 level.")

=== Klasifikasi Stabilitas (2024-01-15 00 UTC) ===
Total grid points       : 12765
Absolutely stable       : 12764  (100.0%)
Conditionally unstable  :     1  (0.0%)
Absolutely unstable     :     0  (0.0%)

Γd digunakan = 9.76 K/km  |  Γm digunakan = 6.0 K/km

=== CAPE Proxy Dua-Level (J/kg) ===
  (Parsel dry-adiabatik dari 850 hPa, dibandingkan environment di 500 hPa)
  Min   : -3568.1
  P5    : -3308.3
  P25   : -3191.8
  Mean  : -3112.7
  Median: -3109.9
  P75   : -3031.6
  P95   : -2917.0
  Max   : -2496.2

  Fraksi kolom dengan CAPE proxy > 0 J/kg : 0.0%
  Catatan: CAPE operasional butuh sounding penuh + LFC/EL, bukan hanya 2 level.

CAPE proxy negatif di semua grid — mean \(\approx -3113\ \text{J/kg}\), maximum pun hanya \(-2496\ \text{J/kg}\) — berarti parsel dry-adiabatik tiba di 500 hPa lebih dingin dari environment, sehingga buoyancy negatif. Ini konsisten secara fisik: parsel mendingin pada \(\Gamma_d \approx 9{,}76\ \text{K/km}\), jauh lebih cepat dari \(\Gamma \approx 5{,}1\ \text{K/km}\) lingkungannya. Setelah naik \({\sim}4{,}4\ \text{km}\), parsel sudah \({\sim}21\ \text{K}\) lebih dingin.

Pelajaran utamanya: updraft cumulonimbus tropis tidak "gratis" dalam kerangka dry-adiabatic. Panas laten kondensasi membuat parsel moist tetap lebih hangat dari environment, sehingga CAPE operasional yang memakai moist adiabat bisa positif meski proxy dry-adiabatic kita selalu negatif.

Diagram di bawah mengilustrasikan mengapa lifting dari 850 hPa bisa menghasilkan buoyancy positif atau negatif tergantung kemiringan \(\Gamma\) relatif terhadap \(\Gamma_d\).

Skema perbandingan suhu parsel (dry adiabat) dan suhu environment di 500 hPa untuk menentukan sign buoyancy.

Hasil Real 2024 dan Langkah Selanjutnya

Alur pipeline tutorial dari cache ke hasil akhir.

Environmental lapse rate 500–850 hPa di atas Indonesia pada 15 Januari 2024 rata-rata \(5{,}11\ \text{K/km}\) — praktis seluruh grid (99,99%) absolutely stable, CAPE proxy \(-3113\ \text{J/kg}\). Studi BMKG untuk Bandung Basin mendokumentasikan CAPE pre-storm \(>1{,}000\ \text{J/kg}\) menjelang hujan lebat — angka yang hanya tercapai dengan kalkulasi moist-adiabatic dan sounding multi-level. Langkah lanjutan:

1. CAPE dan CIN penuh dengan MetPy: gunakan metpy.calc.cape_cin() yang menerima sounding vertikal dengan banyak level. File ERA5 pressure-level dengan 37 level akan menghasilkan CAPE yang lebih akurat daripada pendekatan dua-level kita.

2. Time series harian: terapkan kalkulasi ini ke seluruh 365 hari 2024 untuk satu grid point (misalnya, koordinat Bandung atau Jakarta), lalu korelasikan dengan curah hujan observasi dari BMKG atau IMERG untuk melihat apakah peningkatan \(\Gamma\) mendahului hujan lebat.

3. MALR yang lebih realistis: gunakan data specific humidity dari /data/era5/era5_q_pl500-850_indonesia_2024_d.nc untuk menghitung moist adiabatic lapse rate aktual (\(\Gamma_m\)) secara numerik, menggantikan aproksimasi statis 6 K/km dengan nilai yang bervariasi sesuai suhu dan tekanan.

4. Perbandingan musiman: bandingkan distribusi \(\Gamma\) antara DJF (musim hujan Jawa) dan JJA (musim kering) untuk mengkuantifikasi variabilitas stabilitas atmosfer antar-musim di Indonesia.

Eksplorasi artikel meteorologi lainnya di meteo.my.id — kunjungi meteo.my.id untuk topik cuaca dan iklim Indonesia lainnya.

Referensi

• Lapse Rate — Storm Prediction Center, NOAA — Definisi environmental lapse rate dan tiga kelas stabilitas atmosfer dengan threshold kuantitatif \(\Gamma_d\) dan \(\Gamma_m\); NOAA SPC.
• Explanation of SPC Severe Weather Parameters — NOAA Storm Prediction Center — Definisi CAPE dan CIN, threshold operasional instabilitas, dan hubungan lapse rate dengan potensi konvektif; NOAA SPC.
• ERA5: Data Documentation — ECMWF/Copernicus Knowledge Base — Dokumentasi resmi ERA5 mencakup unit geopotential (\(\text{m}^2/\text{s}^2\)), konversi ke geopotential height (m), dan pressure level yang tersedia; ECMWF, 2024.
• ERA5: Compute Pressure and Geopotential on Model Levels, Geopotential Height and Geometric Height — ECMWF Confluence — Formula presisi konversi geopotential ke geopotential height (\(h = \Phi / 9{,}80665\)) dan perbedaan antara geopotential height dan geometric altitude; ECMWF.
• Study of Atmospheric Stability Indices Forecast for Identifying Convective Activity in the Bandung Basin — Jurnal Meteorologi dan Geofisika (BMKG) — Studi BMKG menerapkan CAPE, K Index, Total Totals, dan Lifted Index dari simulasi WRF-ARW untuk kejadian banjir bandang di Jawa Barat, menunjukkan tren naik indeks stabilitas 2–5 jam sebelum onset; Vol. 26 No. 1, 2026.