论文复现 | 第 04 章：海洋动力学——风应力旋度与海面高度

用户11172986

发布于 2026-04-24 19:35:08

900

第 04 章：海洋动力学——风应力旋度与海面高度

1. 科学背景

ENSO 的演变不仅仅发生在海表面。风场驱动的海洋动力响应才是 La Niña 循环的核心。Figure 6 关注的是风应力旋度 (WSC) 和 海面高度 (SSH) 的演变。

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import warnings
warnings.filterwarnings('ignore')

# --- 配置 ---
DATA_DIR = r"./data"
SSH_FILE = f"{DATA_DIR}/SODA_ssh_mon_1948-2010.nc"
TAUX_FILE = f"{DATA_DIR}/uflx.mon.mean.nc"
TAUY_FILE = f"{DATA_DIR}/vflx.mon.mean.nc"

# 事件列表 (参考表 1)
SINGLE_YEAR_EVENTS = [1964, 1988, 1995, 2005]
MULTI_YEAR_EVENTS = [1949, 1954, 1970, 1973, 1998, 2007]

# 气候态周期与范围
CLIM_START = '1981-01-01'
CLIM_END = '2010-12-31'
LON_RANGE = slice(120, 280)  # 120E to 80W

2. 核心函数：异常值与 WSC 计算

我们需要计算相对于 1981-2010 气候态的月异常。同时，风应力旋度的计算涉及球面坐标系下的偏导数。

def calculate_anomaly(ds, var_name):
    """计算月异常"""
    clim = ds.sel(time=slice(CLIM_START, CLIM_END)).groupby("time.month").mean("time")
    anom = ds.groupby("time.month") - clim
    return anom[var_name]

def calculate_curl(taux, tauy):
    """计算风应力旋度: curl = d(tauy)/dx - d(taux)/dy"""
    R = 6.371e6
    lat, lon = taux.lat.values, taux.lon.values
    lat_rad, lon_rad = np.deg2rad(lat), np.deg2rad(lon)
    
    # 梯度计算
    dtauy_dx = np.gradient(tauy.values, lon_rad, axis=-1) / (R * np.cos(lat_rad)[None, :, None])
    dtaux_dy = np.gradient(taux.values, lat_rad, axis=-2) / R
    
    return xr.DataArray(dtauy_dx - dtaux_dy, coords=taux.coords, dims=taux.dims, name='curl')

3. 合成分析 (Composite Analysis)

我们将不同年份的 La Niña 事件进行合成，展示从发展期 (Apr0) 到衰减期 (Oct1) 的平均演变。

def get_composite_hovmoller(var, events, n_months=19):
    """创建合成 Hovmoller 矩阵 (19个月跨度)"""
    composites = []
    for year in events:
        start_date = f"{year}-04-01"
        end_date = f"{year+1}-10-31"
        try:
            chunk = var.sel(time=slice(start_date, end_date))
            if len(chunk.time) >= 18:
                vals = chunk.values[:19]
                if len(vals) < 19: # 补齐
                    padded = np.full((19, *vals.shape[1:]), np.nan)
                    padded[:len(vals)] = vals
                    vals = padded
                composites.append(vals)
        except:
            continue
    return np.nanmean(composites, axis=0) if composites elseNone