Mengapa Thermal Wind Menautkan Suhu dan Angin
Shear vertikal angin geostrofik bukan hal yang acak — ia ditentukan sepenuhnya oleh gradien suhu horizontal di lapisan atmosfer yang sama. Thermal wind adalah vektor perbedaan antara angin geostrofik di dua level tekanan, dan arah serta besarnya mencerminkan distribusi suhu di antara dua level tersebut. Jet stream subtropis dan polar terbentuk tepat di atas kontur ketebalan yang rapat — daerah di mana gradien suhu horizontal sangat kuat antara tropis yang hangat dan kutub yang dingin.
Sumber: NASA/GSFC Scientific Visualization Studio (The Polar Jet Stream)
Pertanyaan yang kita eksplorasi dalam tutorial ini: dapatkah kita menurunkan shear angin dari pola ketebalan di data ERA5? Dan apa yang terjadi jika kita menerapkan perhitungan ini di lintang rendah seperti Indonesia, di mana asumsi geostrofik tidak berlaku?
Bentuk terintegrasi thermal wind antara dua level tekanan adalah:
$$\mathbf{V}_T = \mathbf{v}_{g,\,p_2} - \mathbf{v}_{g,\,p_1} = \frac{R_d}{f}\,\hat{k} \times \nabla_p \bar{T}\,\ln\!\left(\frac{p_1}{p_2}\right)$$
di mana \(p_1 = 850\ \text{hPa}\) (level bawah), \(p_2 = 500\ \text{hPa}\) (level atas), \(f\) adalah parameter Coriolis, dan \(\nabla_p \bar{T}\) adalah gradien suhu rata-rata lapisan pada permukaan konstan tekanan.
Menyiapkan Data ERA5 di 500 dan 850 hPa
ERA5 menyimpan geopotensial (z, paramId 129) dalam satuan m²/s², bukan meter. Dataset ini tersedia di resolusi horizontal sekitar \(0{,}28125°\) (\({\sim}31\ \text{km}\)) pada 37 level tekanan termasuk 500 dan 850 hPa. Untuk tutorial ini kita butuh tiga variabel: geopotensial (z), angin zonal (u), dan angin meridional (v) — subset Indonesia (6°N–11°S, 95°E–141°E), sepanjang 2024, pukul 00 UTC.
Ketiga variabel diunduh lewat CDS dengan satu loop. Setiap retrieve() digerbang oleh os.path.exists, jadi menjalankan ulang snippet tidak mengunduh dua kali. Antrian CDS biasanya butuh 1–3 menit per variabel pada pemanggilan pertama:
import os, cdsapi, xarray as xr
# Tiga variabel di dua level tekanan; unduh sekali, lewati jika file sudah ada
VARS = {
"geopotential": "era5_z_pl500-850_indonesia_2024_d.nc",
"u_component_of_wind": "era5_u_pl500-850_indonesia_2024_d.nc",
"v_component_of_wind": "era5_v_pl500-850_indonesia_2024_d.nc",
}
for variable, out in VARS.items():
if not os.path.exists(out):
c = cdsapi.Client(quiet=True)
c.retrieve(
"reanalysis-era5-pressure-levels",
{
"product_type": "reanalysis",
"variable": [variable],
"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,
)
zds = xr.open_dataset(VARS["geopotential"])
uds = xr.open_dataset(VARS["u_component_of_wind"])
vds = xr.open_dataset(VARS["v_component_of_wind"])
print("Dimensi z:", dict(zds.sizes))
print("Pressure levels:", zds["pressure_level"].values.tolist())
Dimensi z: {'valid_time': 366, 'pressure_level': 2, 'latitude': 69, 'longitude': 185}
Pressure levels: [850.0, 500.0]
Output di atas mengonfirmasi dua level tekanan (500 dan 850 hPa) dan 366 step waktu harian sepanjang 2024 tersedia di dataset.
Menghitung Ketebalan Lapisan dari Persamaan Hipsometrik
Persamaan hipsometrik menghubungkan ketebalan lapisan dengan suhu virtual rata-rata lapisan:
$$\Delta z = \frac{R_d}{g}\,\bar{T}_v\,\ln\!\left(\frac{p_1}{p_2}\right) \approx 29{,}3\;\bar{T}_v\,\ln\!\left(\frac{850}{500}\right)$$
Lapisan yang lebih hangat mengembang, sehingga ketebalannya lebih besar. Geopotensial ERA5 dalam m²/s² perlu dikonversi ke ketinggian geopotensial (geopotential metre, gpm) dengan:
$$Z_{\text{height}} = \frac{z_{\text{ERA5}}}{g_0}, \quad g_0 = 9{,}80665\ \text{m\,s}^{-2}$$
Kita menggunakan ketebalan 850–500 hPa alih-alih 1000–500 hPa yang umum di midlatitude, karena permukaan 1000 hPa sering berada di bawah atau mendekati permukaan tanah di sebagian wilayah Indonesia.
import numpy as np
# Rata-rata Juli 2024 sebagai contoh representatif
z_jul = zds["z"].sel(valid_time=slice("2024-07-01", "2024-07-31")).mean("valid_time")
# Ketebalan 850–500 hPa dalam geopotential metres (gpm)
thick = (z_jul.sel(pressure_level=500) - z_jul.sel(pressure_level=850)) / 9.80665
print("Ketebalan 850-500 hPa (gpm):")
print(f" min = {float(thick.min()):.1f} gpm")
print(f" mean = {float(thick.mean()):.1f} gpm")
print(f" max = {float(thick.max()):.1f} gpm")
Ketebalan 850-500 hPa (gpm):
min = 4347.5 gpm
mean = 4365.4 gpm
max = 4374.6 gpm
Patokan 5400 gpm yang lazim di midlatitude merujuk pada ketebalan 1000–500 hPa — lapisan yang lebih tebal dan bukan yang kita hitung di sini. Untuk lapisan 850–500 hPa, nilai tipikal di tropis berada di kisaran 4300–4400 gpm, konsisten dengan troposfer tropis yang hangat. Selisih kecil antara nilai min dan max mencerminkan gradien suhu antar sub-region yang nantinya menggerakkan thermal wind.
Rantai konseptual dari gradien suhu ke thermal wind hingga shear geostrofik vertikal.
Menurunkan Vektor Thermal Wind dari Gradien Ketebalan
Dalam koordinat tekanan, thermal wind terintegrasi diekspresikan lewat gradien geopotensial delta-Z antara dua level:
$$u_T = -\frac{1}{f}\frac{\partial(\Delta\Phi)}{\partial y}, \qquad v_T = \frac{1}{f}\frac{\partial(\Delta\Phi)}{\partial x}$$
di mana \(\Delta\Phi = \Phi_{500} - \Phi_{850}\) adalah selisih geopotensial dalam m²/s². Gradien dihitung numeris menggunakan differentiate xarray, lalu dikonversi ke satuan fisik (m/s² per meter jarak) menggunakan faktor metrik.
Satu catatan kritis: parameter Coriolis \(f = 2\Omega\sin\varphi\) mendekati nol di ekuator. Kita mask band \(|\varphi| < 6°\) karena di dekat ekuator \(f\) begitu kecil sehingga faktor \(1/f\) memperkuat noise gradien menjadi vektor berkecepatan ratusan m/s yang jelas tidak fisis. Di luar band itu, thermal wind kembali masuk akal.
# Delta-Phi dalam m²/s² (geopotensial langsung, tanpa dibagi g)
dz = z_jul.sel(pressure_level=500) - z_jul.sel(pressure_level=850)
# Parameter Coriolis sebagai 1D DataArray (xarray auto-broadcast ke 2D)
Omega = 7.2921e-5
f_1d = 2 * Omega * np.sin(np.deg2rad(dz["latitude"]))
# Radius Bumi
a = 6.371e6 # m
# Gradien fisik menggunakan faktor metrik (derajat -> meter)
cos_lat = np.cos(np.deg2rad(dz["latitude"]))
d_dy = dz.differentiate("latitude") * (180.0 / (np.pi * a))
d_dx = dz.differentiate("longitude") * (180.0 / (np.pi * a * cos_lat))
# Komponen thermal wind
uT = -(1.0 / f_1d) * d_dy
vT = (1.0 / f_1d) * d_dx
# Mask band ekuatorial |lat| < 6 derajat: dekat ekuator f begitu kecil
# sehingga 1/f meledakkan noise gradien jadi vektor tak fisis
eq_mask = np.abs(dz["latitude"]) < 6.0
uT = uT.where(~eq_mask)
vT = vT.where(~eq_mask)
speed = np.hypot(uT.values, vT.values)
print(f"Kecepatan thermal wind (m/s):")
print(f" nanmean = {np.nanmean(speed):.2f} m/s")
print(f" nanmax = {np.nanmax(speed):.2f} m/s")
Kecepatan thermal wind (m/s):
nanmean = 3.81 m/s
nanmax = 32.87 m/s
Setelah masking, kecepatan thermal wind turun ke rentang fisis — rata-rata beberapa m/s dengan maksimum sekitar 30 m/s. Nilai terkuat muncul di pita luar-ekuatorial (mendekati 6°N dan pita 6–11°S), di mana \(f\) sudah cukup besar untuk memberi estimasi yang bermakna sekaligus mencerminkan gradien suhu yang lebih terstruktur.
Membandingkan Thermal Wind dengan Shear Angin Aktual
Shear angin aktual dihitung sebagai selisih komponen angin antara 500 dan 850 hPa. Di midlatitude, thermal wind dan shear geostrofik aktual hampir identik. Di pita luar-ekuatorial Indonesia kita bisa cek seberapa dekat keduanya lewat korelasi Pearson dan RMS, memakai mask \(|\varphi| < 6°\) yang sama seperti sebelumnya.
# Shear angin aktual (rata-rata Juli 2024)
u_jul = uds["u"].sel(valid_time=slice("2024-07-01", "2024-07-31")).mean("valid_time")
v_jul = vds["v"].sel(valid_time=slice("2024-07-01", "2024-07-31")).mean("valid_time")
uS = u_jul.sel(pressure_level=500) - u_jul.sel(pressure_level=850)
vS = v_jul.sel(pressure_level=500) - v_jul.sel(pressure_level=850)
# Terapkan mask ekuatorial yang sama
uS = uS.where(~eq_mask)
vS = vS.where(~eq_mask)
# RMS perbedaan (menggunakan array values agar shape konsisten)
du = uT.values - uS.values
dv = vT.values - vS.values
valid = ~np.isnan(du) & ~np.isnan(dv)
rms = np.sqrt(np.mean(du[valid]**2 + dv[valid]**2))
r_u = np.corrcoef(uT.values[valid], uS.values[valid])[0, 1]
r_v = np.corrcoef(vT.values[valid], vS.values[valid])[0, 1]
print(f"RMS perbedaan thermal wind vs shear aktual: {rms:.2f} m/s")
print(f"Korelasi Pearson komponen u: {r_u:.3f}")
print(f"Korelasi Pearson komponen v: {r_v:.3f}")
RMS perbedaan thermal wind vs shear aktual: 4.34 m/s
Korelasi Pearson komponen u: 0.595
Korelasi Pearson komponen v: 0.245
Komponen zonal (\(u\)) berkorelasi cukup kuat — sekitar \(0{,}6\) — dengan shear aktual, menandakan keseimbangan thermal wind sebagian besar berlaku di pita ini. Komponen meridional (\(v\)) jauh lebih lemah (sekitar \(0{,}25\)), mencerminkan pengaruh ageostrofik dan monsoon yang lebih besar plus \(f\) yang masih marginal di lintang rendah. RMS perbedaan sekitar \(4\ \text{m/s}\) — kecil dibanding blow-up sebelum masking, tetapi mengingatkan bahwa thermal wind di tropis adalah aproksimasi, bukan identitas.
Visualisasi Ketebalan dan Thermal Wind di Atas Indonesia
Peta berikut menampilkan ketebalan 850–500 hPa sebagai filled contour dan panah thermal wind sebagai overlay. Pita \(|\varphi| < 6°\) ditutup mask sehingga panah hanya muncul di sabuk luar-ekuatorial tempat thermal wind bermakna.
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
fig = plt.figure(figsize=(11, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([95, 141, -11, 6], crs=ccrs.PlateCarree())
lon2d, lat2d = np.meshgrid(thick["longitude"].values, thick["latitude"].values)
# Filled contour ketebalan
cf = ax.contourf(
lon2d, lat2d, thick.values,
levels=20, cmap="RdYlBu_r",
transform=ccrs.PlateCarree()
)
cb = plt.colorbar(cf, ax=ax, orientation="horizontal", pad=0.05, shrink=0.8)
cb.set_label("Ketebalan 850-500 hPa (gpm)", fontsize=10)
# Thermal wind vectors (subsample setiap 6 titik grid)
step = 6
ax.quiver(
lon2d[::step, ::step],
lat2d[::step, ::step],
uT.values[::step, ::step],
vT.values[::step, ::step],
scale=300, width=0.003, color="black", alpha=0.7,
transform=ccrs.PlateCarree()
)
ax.coastlines("50m", linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle="--")
gl = ax.gridlines(draw_labels=True, linewidth=0.5, color="gray", alpha=0.5, linestyle="--")
gl.top_labels = False
gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(range(95, 145, 10))
gl.ylocator = mticker.FixedLocator(range(-10, 10, 5))
ax.set_title(
"Ketebalan 850-500 hPa dan Thermal Wind — Indonesia, Juli 2024",
fontsize=12
)
plt.tight_layout()
plt.savefig("thermal_wind_indonesia.png", dpi=140, bbox_inches="tight")
print("saved figure")
Warna merah-kuning menunjukkan lapisan lebih tebal (lebih hangat), warna biru lapisan lebih tipis (lebih dingin). Panah hanya muncul di pita luar-ekuatorial (di luar \(|\varphi| < 6°\)) karena di sanalah thermal wind bermakna; panah terpanjang berada di sub-domain dengan gradien ketebalan terkuat.
Kaveat Lintang Rendah dan Langkah Selanjutnya
Parameter Coriolis \(f = 2\Omega\sin\varphi\) mendekati nol tepat di ekuator (\(\varphi = 0°\)). Seluruh konsep thermal wind bertumpu pada keseimbangan geostrofik, yang membutuhkan \(f\) cukup besar agar gaya Coriolis menyeimbangkan gaya tekanan. Di wilayah Indonesia (6°N–11°S) nilai \(|f|\) hanya sekitar \(10^{-5}\ \text{s}^{-1}\) — satu orde lebih kecil dari nilai midlatitude (\(|f| \sim 10^{-4}\ \text{s}^{-1}\) di 45°).
Itulah mengapa kita mask band \(|\varphi| < 6°\) di snippet-3: di dekat ekuator faktor \(1/f\) meledak dan formula thermal wind kehilangan makna fisis. Di luar band itu — di 6°N maupun pita 6–11°S — korelasi komponen zonal yang mencapai \(\sim 0{,}6\) menunjukkan thermal wind sudah cukup andal sebagai proksi gradien suhu, walau komponen meridional tetap lemah. Untuk analisis forecasting di Indonesia, shear angin aktual dari snippet-4 (u500 − u850, v500 − v850) tetap menjadi ukuran yang paling relevan secara operasional.
Beberapa langkah eksplorasi selanjutnya yang menarik: analisis variabilitas ketebalan sepanjang musim 2024 untuk mendeteksi sinyal monsoon, komputasi thermal wind pada domain yang mencakup lintang midlatitude (Australia selatan, 30–40°S) di mana geostrofi lebih valid, atau analisis korelasi antara pola ketebalan dan curah hujan bulanan.
Eksplorasi artikel meteorologi lainnya di meteo.my.id (https://meteo.my.id).
Referensi
- Thermal Wind Effect (Practical Meteorology, Stull) — Definisi dan komponen thermal wind, kaitan dengan veering/backing angin, dan konteks meteorologi sinoptik midlatitude.
- Hypsometric Equation (Practical Meteorology, Stull) — Penurunan persamaan hipsometrik dari hukum gas ideal dan hidrostatis; ketebalan proporsional dengan suhu virtual rata-rata lapisan.
- Constant Pressure Charts: Thickness — NOAA JetStream — Penggunaan peta ketebalan dalam prakiraan sinoptik: nilai 540 dam sebagai batas hujan/salju dan packing kontur sebagai indikator frontal.
- ERA5: compute pressure and geopotential on model levels — Cara mengonversi geopotensial ERA5 (m²/s²) ke ketinggian geopotensial menggunakan \(g_0 = 9{,}80665\ \text{m\,s}^{-2}\).
- ERA5: data documentation — Resolusi spasial ERA5, 37 level tekanan, dan metadata variabel geopotensial (shortName z, paramId 129).
Tidak ada komentar:
Posting Komentar