Menghitung Frekuensi Brunt-Väisälä dari ERA5

Mengapa Frekuensi Brunt-Väisälä Penting di Tropis

Setiap kali pilot melaporkan turbulence di langit yang tampak bersih di atas Kalimantan, atau kita melihat awan lentikular mengambang di puncak gunung Jaya Wijaya, ada satu parameter fisika yang bekerja di baliknya: frekuensi Brunt-Väisälä, \(N\). Parameter ini mengukur seberapa kuat atmosfer "mengembalikan" parsel udara yang terusik ke posisi semula. Semakin besar \(N\), semakin stabil kolom udara dan semakin cepat parsel berosilasi sebelum berhenti.

Di Indonesia, static stability terasa sangat relevan. Maritime continent kita adalah salah satu pusat deep convection terkuat di Bumi — sebagian besar karena stabilitas statik troposfer yang lemah di atas lautan hangat, di mana \(N^2\) kerap mendekati nol. Sebaliknya, lapisan dengan \(N\) besar memungkinkan gravity wave merambat vertikal — kunci memahami gelombang orografis di Papua, dispersi polutan kota, maupun turbulence di jalur penerbangan.

Secara formal, frekuensi Brunt-Väisälä kuadrat untuk udara kering didefinisikan sebagai:

$$N^2 = \frac{g}{\theta} \cdot \frac{d\theta}{dz}$$

dengan \(g = 9{,}80665\ \text{m s}^{-2}\), \(\theta\) adalah potential temperature (K), dan \(z\) adalah ketinggian (m). Jika \(N^2 > 0\), atmosfer stabil dan parsel berosilasi; jika \(N^2 < 0\), parsel berakselerasi menjauhi ekuilibriumnya dan konveksi berkembang.

Dalam tutorial ini kita download data ERA5 pressure-level (temperature dan geopotential di 500 dan 850 hPa) untuk Indonesia tahun 2024, lalu menghitung \(N^2\) bulk lapisan 850–500 hPa langsung dari definisinya dengan Python dan xarray.

Mengunduh dan Memuat Data ERA5

ERA5 menyimpan geopotential (m² s⁻²) dan temperature (K) pada 37 level tekanan standar, resolusi ~31 km, dari 1940 hingga nyaris real-time. Untuk tutorial ini kita butuh dua variabel di 500 dan 850 hPa: temperature dan geopotential.

Sebelum mulai, daftar akun dan dapatkan API key di cds.climate.copernicus.eu, lalu pip install cdsapi xarray netCDF4. Snippet berikut men-download dua file NetCDF untuk Indonesia 2024 hanya satu kali — guard if not os.path.exists melewati download pada run berikutnya. Tergantung antrean CDS, tiap file selesai dalam 1–10 menit.

import os, cdsapi, xarray as xr, numpy as np

OUT_T = "era5_t_pl500-850_indonesia_2024_d.nc"
OUT_Z = "era5_z_pl500-850_indonesia_2024_d.nc"

if not os.path.exists(OUT_T):
    c = cdsapi.Client(quiet=True)
    c.retrieve(
        "reanalysis-era5-pressure-levels",
        {
            "product_type": "reanalysis",
            "variable": ["temperature"],
            "pressure_level": ["500", "850"],
            "year": "2024",
            "month": [f"{m:02d}" for m in range(1, 13)],
            "day":   [f"{d:02d}" for d in range(1, 32)],
            "time":  "00:00",
            "area":  [6, 95, -11, 141],
            "format": "netcdf",
        },
        OUT_T,
    )

if not os.path.exists(OUT_Z):
    c = cdsapi.Client(quiet=True)
    c.retrieve(
        "reanalysis-era5-pressure-levels",
        {
            "product_type": "reanalysis",
            "variable": ["geopotential"],
            "pressure_level": ["500", "850"],
            "year": "2024",
            "month": [f"{m:02d}" for m in range(1, 13)],
            "day":   [f"{d:02d}" for d in range(1, 32)],
            "time":  "00:00",
            "area":  [6, 95, -11, 141],
            "format": "netcdf",
        },
        OUT_Z,
    )

ds_t = xr.open_dataset(OUT_T)
ds_z = xr.open_dataset(OUT_Z)

print("=== Temperature dataset ===")
print(ds_t.sizes)
print("Data vars:", list(ds_t.data_vars))
print("Pressure levels:", ds_t["pressure_level"].values)

print("\n=== Geopotential dataset ===")
print(ds_z.sizes)
print("Data vars:", list(ds_z.data_vars))
print("Pressure levels:", ds_z["pressure_level"].values)

=== Temperature dataset ===
Frozen({'valid_time': 366, 'pressure_level': 2, 'latitude': 69, 'longitude': 185})
Data vars: ['t']
Pressure levels: [850. 500.]

=== Geopotential dataset ===
Frozen({'valid_time': 366, 'pressure_level': 2, 'latitude': 69, 'longitude': 185})
Data vars: ['z']
Pressure levels: [850. 500.]

Output di atas menampilkan dimensi kedua dataset: valid_time (366 langkah harian, 2024 tahun kabisat), pressure_level ([850.0, 500.0]), serta grid 69×185 titik. Temperature tersimpan sebagai t, geopotential sebagai z — keduanya jadi input semua snippet berikutnya.

Menghitung Suhu Potensial di Setiap Level

Potential temperature \(\theta\) adalah suhu yang akan dimiliki parsel jika diturunkan secara adiabatik ke level referensi 1000 hPa. Berbeda dari suhu aktual \(T\) yang turun terhadap ketinggian, \(\theta\) justru naik dengan ketinggian di atmosfer stabil — itulah mengapa gradient \(d\theta/dz > 0\) menjadi kriteria stabilitas statik.

Formulanya untuk udara kering:

$$\theta = T \cdot \left(\frac{1000}{p}\right)^{0{,}286}$$

dengan \(p\) dalam hPa dan \(T\) dalam K. Eksponen \(0{,}286 \approx R_a / c_p\) menggunakan \(R_a = 287\ \text{J kg}^{-1}\text{K}^{-1}\) dan \(c_p = 1005\ \text{J kg}^{-1}\text{K}^{-1}\) — nilai konvensi WMO yang dipakai secara luas.

Snippet berikut menghitung \(\theta\) di 500 dan 850 hPa menggunakan variabel t dari dataset, lalu mencetak domain-mean dan selisihnya.

import numpy as np

# Potential temperature: theta = T * (1000/p)^0.286
theta500 = ds_t["t"].sel(pressure_level=500) * (1000 / 500) ** 0.286
theta850 = ds_t["t"].sel(pressure_level=850) * (1000 / 850) ** 0.286

# Annual-mean domain-mean (average over latitude, longitude, dan valid_time)
t500_mean = float(theta500.mean(["latitude", "longitude", "valid_time"]))
t850_mean = float(theta850.mean(["latitude", "longitude", "valid_time"]))

print(f"Domain-mean θ @ 500 hPa : {t500_mean:.2f} K")
print(f"Domain-mean θ @ 850 hPa : {t850_mean:.2f} K")
print(f"Δθ (500 − 850 hPa)       : {t500_mean - t850_mean:.2f} K")

Domain-mean θ @ 500 hPa : 327.87 K
Domain-mean θ @ 850 hPa : 305.38 K
Δθ (500 − 850 hPa)       : 22.48 K

Output menampilkan dua nilai \(\theta\) dan selisihnya \(\Delta\theta\). Yang penting: \(\theta_{500}\) lebih tinggi dari \(\theta_{850}\), artinya potential temperature naik seiring ketinggian — tanda klasik stabilitas statik. Selisih \(\Delta\theta\) positif inilah pembilang dalam formula \(N^2\) berikutnya.

Frekuensi Brunt-Väisälä dan Periode Osilasi Buoyancy

Untuk menghitung \(N^2\) bulk antara 850 dan 500 hPa, kita gunakan aproksimasi finite-difference dari definisi kontinu \(N^2 = (g/\theta)(d\theta/dz)\):

$$N^2 = \frac{g}{\bar{\theta}} \cdot \frac{\Delta\theta}{\Delta H}$$

dengan \(\bar{\theta} = (\theta_{500} + \theta_{850})/2\) adalah rata-rata layer, \(\Delta\theta = \theta_{500} - \theta_{850}\) (K), dan \(\Delta H = H_{500} - H_{850}\) (m) adalah selisih geopotential height. ERA5 menyimpan geopotential \(\Phi\) dalam satuan m² s⁻², sehingga konversi ke ketinggian dilakukan dengan \(H = \Phi / g_0\), dengan \(g_0 = 9{,}80665\ \text{m s}^{-2}\).

Dari \(N^2\), kita dapatkan \(N = \sqrt{N^2}\) (s⁻¹) dan periode osilasi buoyancy:

$$T_{BV} = \frac{2\pi}{N}\ \text{(detik)}$$

Nilai referensi troposfer bebas yang umum dikutip adalah \(N \approx 0{,}01\ \text{s}^{-1}\), \(N^2 \approx 1\times10^{-4}\ \text{s}^{-2}\), dan \(T_{BV} \approx 10\ \text{menit}\). Di stratosfer yang lebih stabil, \(N\) bisa mencapai \(\approx 0{,}02\ \text{s}^{-1}\) dengan \(T_{BV} \approx 5\ \text{menit}\).

import numpy as np

g = 9.80665  # m/s²

# Convert geopotential Phi (m² s⁻²) to geopotential height H (m)
H500 = ds_z["z"].sel(pressure_level=500) / g
H850 = ds_z["z"].sel(pressure_level=850) / g

# Potential temperature di setiap level
theta500 = ds_t["t"].sel(pressure_level=500) * (1000 / 500) ** 0.286
theta850 = ds_t["t"].sel(pressure_level=850) * (1000 / 850) ** 0.286

# Bulk layer quantities
dtheta     = theta500 - theta850           # K  — positif → stabil
dH         = H500 - H850                   # m  — positif (500 hPa lebih tinggi)
theta_mean = (theta500 + theta850) / 2     # K

# N² = (g / theta_mean) * (dtheta / dH)
N2 = (g / theta_mean) * (dtheta / dH)

# Domain-mean dan annual-mean
N2_dmean = float(N2.mean(["latitude", "longitude", "valid_time"]))
N_dmean  = float(np.sqrt(max(N2_dmean, 0.0)))
T_BV_min = (2 * np.pi / N_dmean) / 60    # detik → menit

print(f"Domain-mean N²   : {N2_dmean:.3e} s⁻²")
print(f"Domain-mean N    : {N_dmean:.5f} s⁻¹")
print(f"Buoyancy period  : {T_BV_min:.1f} menit")

Domain-mean N²   : 1.593e-04 s⁻²
Domain-mean N    : 0.01262 s⁻¹
Buoyancy period  : 8.3 menit

Output di atas menampilkan \(N^2 \approx 1{,}6\times10^{-4}\ \text{s}^{-2}\), \(N \approx 0{,}013\ \text{s}^{-1}\), dan \(T_{BV} \approx 8\ \text{menit}\) untuk lapisan 850–500 hPa Indonesia. Nilai ini sedikit di atas referensi troposfer bebas (\(N^2 \approx 1\times10^{-4}\ \text{s}^{-2}\), \(T_{BV} \approx 10\ \text{menit}\)) — wajar, karena rata-rata bulk 850–500 hPa turut memuat mid-troposfer yang lebih stabil. Periode buoyancy sekitar 8 menit berarti parsel yang terusik berosilasi naik-turun kira-kira sekali tiap 8 menit sebelum teredam.

Ketiga rezim yang ditentukan oleh tanda \(N^2\) dapat dirangkum sebagai berikut:

Diagram 1. Tiga rezim stabilitas atmosfer berdasarkan tanda N². Nilai positif memungkinkan gravity waves dan awan lentikular; nilai negatif memicu deep convection.

Visualisasi N² di Seluruh Indonesia

Angka skalar domain-mean memberi gambaran rata-rata, tetapi pola spasial dan variasi temporal \(N^2\) jauh lebih informatif. Kita buat dua visualisasi yang saling melengkapi: peta annual-mean di seluruh Indonesia dan time series harian sepanjang 2024.

Snippet berikut menghasilkan peta \(N^2\) dengan cartopy. Pola spasialnya memperlihatkan di mana stabilitas troposfer paling kuat — biasanya di atas perairan ber-SST lebih rendah atau wilayah yang dipengaruhi aliran udara benua yang kering.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

g = 9.80665

theta500 = ds_t["t"].sel(pressure_level=500) * (1000 / 500) ** 0.286
theta850 = ds_t["t"].sel(pressure_level=850) * (1000 / 850) ** 0.286
H500 = ds_z["z"].sel(pressure_level=500) / g
H850 = ds_z["z"].sel(pressure_level=850) / g

dtheta     = theta500 - theta850
dH         = H500 - H850
theta_mean = (theta500 + theta850) / 2
N2         = (g / theta_mean) * (dtheta / dH)

# Annual-mean spatial field
N2_annual = N2.mean("valid_time")

fig, ax = plt.subplots(
    1, 1, figsize=(11, 6),
    subplot_kw={"projection": ccrs.PlateCarree()}
)
ax.set_extent([95, 141, -11, 6], crs=ccrs.PlateCarree())

pcm = ax.pcolormesh(
    N2_annual["longitude"].values,
    N2_annual["latitude"].values,
    N2_annual.values,
    transform=ccrs.PlateCarree(),
    cmap="viridis",
)
ax.add_feature(cfeature.COASTLINE.with_scale("50m"), linewidth=0.8, color="white")
gl = ax.gridlines(draw_labels=True, linewidth=0.5, color="gray", alpha=0.5)
gl.top_labels = False
gl.right_labels = False

cbar = fig.colorbar(pcm, ax=ax, orientation="horizontal", pad=0.08, shrink=0.8)
cbar.set_label("N² (s⁻²)", fontsize=12)

ax.set_title(
    "Frekuensi Brunt-Väisälä² — Rata-rata Tahunan 850–500 hPa (ERA5 2024)",
    fontsize=12
)
plt.savefig("bv_map.png", dpi=150, bbox_inches="tight")
print("Saved bv_map.png")

Peta di atas menampilkan variasi spasial \(N^2\) tahunan di atas Indonesia. Warna lebih terang menandai stabilitas troposfer lebih kuat; warna gelap menandai lapisan lebih netral yang cenderung konduktif untuk deep convection.

Snippet kedua memperlihatkan bagaimana \(N^2\) domain-mean berfluktuasi dari hari ke hari sepanjang 2024 — sinyal musiman maupun variabilitas synoptik tercermin di dalamnya.

import matplotlib.pyplot as plt
import numpy as np

g = 9.80665

theta500 = ds_t["t"].sel(pressure_level=500) * (1000 / 500) ** 0.286
theta850 = ds_t["t"].sel(pressure_level=850) * (1000 / 850) ** 0.286
H500 = ds_z["z"].sel(pressure_level=500) / g
H850 = ds_z["z"].sel(pressure_level=850) / g

dtheta     = theta500 - theta850
dH         = H500 - H850
theta_mean = (theta500 + theta850) / 2
N2         = (g / theta_mean) * (dtheta / dH)

# Daily domain-mean (rata-rata spasial atas latitude dan longitude)
N2_ts = N2.mean(["latitude", "longitude"])

fig, ax = plt.subplots(figsize=(11, 4))
ax.plot(
    N2_ts["valid_time"].values,
    N2_ts.values,
    color="#1E5A8A",
    linewidth=0.9,
    label="N² domain-mean harian"
)
ax.axhline(
    1e-4,
    color="orange",
    linestyle="--",
    linewidth=1.2,
    label="Referensi troposfer bebas N² = 1×10⁻⁴ s⁻²"
)
ax.set_xlabel("Tanggal (2024)", fontsize=11)
ax.set_ylabel("N² (s⁻²)", fontsize=11)
ax.set_title("Domain-mean N² Harian — Indonesia 850–500 hPa (ERA5 2024)", fontsize=12)
ax.legend(fontsize=10)
ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig("bv_timeseries.png", dpi=150, bbox_inches="tight")
print("Saved bv_timeseries.png")

Time series di atas memperlihatkan fluktuasi stabilitas troposfer Indonesia sepanjang tahun — perhatikan sinyal musiman antara musim hujan (Oktober–April) dan kemarau (Mei–September) maupun variasi synoptik berperiode pendek. Garis referensi oranye di \(N^2 = 1\times10^{-4}\ \text{s}^{-2}\) mempermudah perbandingan dengan troposfer ekstratropis.

Interpretasi dalam Konteks Meteorologi Tropis

Angka \(N^2\) ini perlu dibaca dalam konteks troposfer tropis. Di ekstratropis, static stability umumnya lebih tinggi karena lapse rate teramati jauh di bawah dry adiabatic lapse rate (\(\Gamma_d \approx 9{,}8\ \text{K/km}\)), sehingga \(d\theta/dz\) besar. Di atas Indonesia — terutama di lapisan bawah dekat permukaan — lapse rate jauh lebih dekat ke moist adiabatic (~6–7 K/km), sehingga \(d\theta/dz\) mengecil dan \(N^2\) turun, kadang mendekati nol saat konveksi aktif.

Inilah fondasi termodinamika mengapa deep convection begitu umum di atas Sumatera, Kalimantan, dan Papua: dengan troposfer yang hampir netral, parsel panas-lembap dari laut hangat (SST > 28°C) hanya perlu sedikit dorongan untuk melampaui level kondensasi dan terus naik tanpa hambatan buoyancy berarti. Di atas warm pool Samudra Pasifik barat, \(N^2\) kerap mendekati nol atau sesaat negatif — pemicu awal badai konvektif skala meso.

Skenario kontras berlaku di area dengan static stability lebih tinggi. Saat aliran udara stabil dipaksa naik oleh orografi — misalnya melintasi pegunungan Jaya Wijaya (>4.000 m) — nilai \(N\) yang signifikan melahirkan gravity wave yang efisien. Gelombang ini merambat ke atas dan di permukaan tampak sebagai awan lentikular di sisi lee gunung atau gelombang orografis pada citra satellite inframerah.

Ada hubungan mendasar di sini: tidak ada internal gravity wave yang berosilasi lebih cepat dari \(N\) lokal. Lapisan dengan \(N\) besar memuat rentang periode gravity wave yang luas, sedangkan lapisan mendekati netral hanya meloloskan gelombang berperiode sangat panjang. Memetakan distribusi \(N^2\) dari ERA5 karena itu menjadi langkah pertama menilai di mana gravity wave dapat berkembang di atas Indonesia.

Langkah Selanjutnya

Tutorial ini menghitung \(N^2\) sebagai aproksimasi bulk dua level (850–500 hPa). Ada beberapa arah eksplorasi yang natural dari sini:

1. Profil vertikal \(N^2\): gunakan semua 37 level tekanan ERA5 untuk mendapatkan \(N^2(p)\) sebagai fungsi ketinggian — dari boundary layer hingga lower stratosphere — dan lihat lapisan tropopause menonjol sebagai lompatan stabilitas tajam.
2. Perbandingan troposfer vs stratosfer: tambahkan level 50–200 hPa dan konfirmasi \(N^2\) stratosfer jauh lebih tinggi daripada troposfer.
3. Scorer parameter dan wave trapping: kombinasikan \(N\) dengan profil shear angin untuk menghitung Scorer parameter \(l^2 = N^2/U^2 - (d^2U/dz^2)/U\) dan nilai kondisi propagasi vs trapping untuk gravity wave.
4. Validasi dengan radiosonde: bandingkan \(N^2\) ERA5 dengan data radiosonde BMKG (misalnya Koto Tabang atau Biak) untuk mengukur seberapa akurat reanalysis mereproduksi stabilitas troposfer aktual.

Frekuensi Brunt-Väisälä adalah parameter fundamental dinamika atmosfer. Begitu bisa menghitungnya dari ERA5, analisis gravity wave, parametrisasi konveksi, dan kajian turbulence pun terbuka lebar.

Eksplorasi artikel meteorologi lainnya di meteo.my.id — kunjungi https://meteo.my.id untuk artikel lain seputar cuaca dan data atmosfer.

Referensi

• ERA5: data documentation – Copernicus Knowledge Base, ECMWF — Dokumentasi resmi ERA5 yang mencakup resolusi temporal dan spasial (~31 km), 37 level tekanan standar, serta definisi variabel pressure-level termasuk geopotential (m² s⁻²) dan temperature (K).
• ERA5: compute geopotential height and geometric height – ECMWF Confluence Wiki — Halaman teknis ECMWF yang mengkonfirmasi konversi geopotential ke geopotential height menggunakan konstanta gravitasi standar \(g_0 = 9{,}80665\ \text{m s}^{-2}\).
• The Brunt-Väisälä Frequency – Gravity Waves (EUMETRAIN / COMET Program) — Sumber pelatihan EUMETSAT yang mendefinisikan \(N^2\) dalam konteks atmosferik dan menghubungkannya dengan mountain wave, gravity wave, serta awan lentikular sebagai manifestasi fisik stabilitas.
• Brunt–Väisälä frequency – Wikipedia — Artikel ensiklopedia yang merangkum formulasi \(N^2\) dalam potential temperature maupun suhu aktual, beserta interpretasi tiga rezim stabilitas (\(N^2 > 0\), \(= 0\), \(< 0\)) dan hubungannya dengan internal gravity wave.
• Buoyancy (gravity) waves in the atmosphere – University of Western Ontario, Dr. W. Hocking — Catatan kuliah yang memberikan formula \(T_{BV} = 2\pi/N\) dengan nilai tipikal troposfer (~10 menit, \(N \approx 0{,}01\ \text{s}^{-1}\)) dan stratosfer (~5 menit, \(N \approx 0{,}02\ \text{s}^{-1}\)).