Calculate the SM-T coupling using ifs_tco3999-ng5_rcbmf_cf¶

Author: Yifan Cheng

Reference:

Miralles et al. (2012)

import datetime
current_time = datetime.datetime.now().strftime("%Y-%m-%d %H:%M")
print(f'Last updated at {current_time}')

Last updated at 2025-05-16 09:43

from pathlib import Path 
import xarray as xr
import cartopy.crs as ccrs
import uxarray as ux
import numpy as np
import cartopy.feature as cf
import matplotlib.pyplot as plt
import pandas as pd

import intake

cat_url = "https://digital-earths-global-hackathon.github.io/catalog/catalog.yaml"
cat = intake.open_catalog(cat_url).NCAR

# check on global models
model_run2 = cat['ifs_tco3999-ng5_rcbmf_cf']

ds_fine2 = model_run2(zoom=7,time='PT1H').to_dask()

/glade/work/yifanc17/miniconda3/envs/hackathon_env/lib/python3.12/site-packages/intake_xarray/base.py:21: FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension names in future, in order to be more consistent with `DataArray.dims`. To access a mapping from dimension names to lengths, please use `Dataset.sizes`.
  'dims': dict(self._ds.dims),

uxds_fine2 = ux.UxDataset.from_healpix(ds_fine2)

# for var_name, da in uxds_fine2.data_vars.items():
#     if da.dims == ('time', 'n_face') and da.isel(time=0).isnull().all():
#         print(f"Variable '{var_name}' is all NaN")

var_longname_dict2 = {
    var: ds_fine2[var].attrs.get("name", "N/A")
    for var in ds_fine2.data_vars
}

var_longname_dict2

{'100si': '100 metre wind speed',
 '100u': '100 metre U wind component',
 '100v': '100 metre V wind component',
 '10si': '10 metre wind speed',
 '2d': '2 metre dewpoint temperature',
 'avg_sd24': 'Time-mean snow depth water equivalent',
 'blh': 'Boundary layer height',
 'cc': 'Fraction of cloud cover',
 'chnk': 'Charnock',
 'ciwc': 'Specific cloud ice water content',
 'clivi': 'Total column cloud ice water',
 'clt': 'Total cloud cover',
 'clwc': 'Specific cloud liquid water content',
 'clwvi': 'Total column cloud liquid water',
 'crwc': 'Specific rain water content',
 'cswc': 'Specific snow water content',
 'e': 'Evaporation',
 'fdir': 'Total sky direct short-wave (solar) radiation at surface',
 'hcc': 'High cloud cover',
 'hflsd': 'Surface latent heat flux',
 'hfssd': 'Surface sensible heat flux',
 'hur': 'Relative humidity',
 'hus': 'Specific humidity',
 'i10fg': 'Instantaneous 10 metre wind gust',
 'lcc': 'Low cloud cover',
 'lgws': 'Eastward gravity wave surface stress',
 'litota1': 'Averaged total lightning flash density in the last hour',
 'litoti': 'Instantaneous total lightning flash density',
 'lsm': 'Land-sea mask',
 'lsp': 'Large-scale precipitation',
 'lwcs': 'Liquid water content in snow pack',
 'mcc': 'Medium cloud cover',
 'mgws': 'Northward gravity wave surface stress',
 'mrso': 'Volumetric soil water layer 1',
 'mtpr': 'Mean total precipitation rate',
 'mucape': 'Most-unstable CAPE',
 'pr': 'Total precipitation',
 'prs': 'Snowfall',
 'prw': 'Total column vertically-integrated water vapour',
 'ps': 'Surface pressure',
 'psl': 'Mean sea level pressure',
 'pv': 'Potential vorticity',
 'qall': 'Specific rain water content',
 'rlds': 'Surface long-wave (thermal) radiation downwards',
 'rlus': 'Surface long-wave (thermal) radiation downwards',
 'rluscs': 'Surface net long-wave (thermal) radiation, clear sky',
 'rlut': 'Top net long-wave (thermal) radiation',
 'rlutcs': 'Top net long-wave (thermal) radiation, clear sky',
 'rsds': 'Surface short-wave (solar) radiation downwards',
 'rsdscs': 'Surface short-wave (solar) radiation downward clear-sky',
 'rsdt': 'TOA incident short-wave (solar) radiation',
 'rsn': 'Snow density',
 'rsus': 'Surface short-wave (solar) radiation downwards',
 'rsuscs': 'Surface short-wave (solar) radiation downward clear-sky',
 'rsut': 'TOA incident short-wave (solar) radiation',
 'sd': 'Snow depth water equivalent',
 'siconc': 'Sea ice area fraction',
 'sro': 'Surface runoff',
 'ssro': 'Sub-surface runoff',
 'sst': 'Sea surface temperature',
 'stl1': 'Soil temperature level 1',
 'stl2': 'Soil temperature level 2',
 'stl3': 'Soil temperature level 3',
 'stl4': 'Soil temperature level 4',
 'swe': 'Snow depth',
 'ta': 'Temperature',
 'tas': '2 metre temperature',
 'tauu': 'Time-integrated eastward turbulent surface stress',
 'tauv': 'Time-integrated northward turbulent surface stress',
 'tcrw': 'Total column rain water',
 'tcsw': 'Total column snow water',
 'tprate': 'Total precipitation rate',
 'ts': 'Skin temperature',
 'tsrc': 'Top net short-wave (solar) radiation, clear sky',
 'ua': 'U component of wind',
 'uas': '10 metre U wind component',
 'va': 'V component of wind',
 'vas': '10 metre V wind component',
 'wa': 'Vertical velocity',
 'z': 'Geopotential',
 'zg': 'Geopotential'}

#slhf #Surface latent heat flux
#sshf #Sensible heat flux
#ssr #Surface net short-wave (solar) radiation
#str #Surface net long-wave (thermal) radiation
#2t #2 metre temperatur

land_mask=ux.UxDataArray((uxds_fine2.sst.isel(time=0)==9999).astype(int),uxgrid=uxds_fine2.uxgrid)

sel_vswc=(uxds_fine2['mrso']).where(land_mask)
sel_le=(uxds_fine2['hflsd']).where(land_mask) #W/m2
sel_h=(uxds_fine2['hfssd']).where(land_mask) #W/m2
sel_rnet=((uxds_fine2['rlds']-uxds_fine2['rlus']+uxds_fine2['rsds']-uxds_fine2['rsus'])).where(land_mask) #W/m2
sel_t2m=(uxds_fine2['tas']-273.15).where(land_mask) #C
sel_sp=(uxds_fine2['ps']/1000).where(land_mask) #kPa

sel_vswc_DJF=sel_vswc.sel(time=(uxds_fine2.time.dt.month.isin([1,2,12]))&(uxds_fine2.time.dt.year == 2020))
sel_le_DJF=sel_le.sel(time=(uxds_fine2.time.dt.month.isin([1,2,12]))&(uxds_fine2.time.dt.year == 2020))
sel_h_DJF=sel_h.sel(time=(uxds_fine2.time.dt.month.isin([1,2,12]))&(uxds_fine2.time.dt.year == 2020))
sel_rnet_DJF=sel_rnet.sel(time=(uxds_fine2.time.dt.month.isin([1,2,12]))&(uxds_fine2.time.dt.year == 2020))
sel_t2m_DJF=sel_t2m.sel(time=(uxds_fine2.time.dt.month.isin([1,2,12]))&(uxds_fine2.time.dt.year == 2020))
sel_sp_DJF=sel_sp.sel(time=(uxds_fine2.time.dt.month.isin([1,2,12]))&(uxds_fine2.time.dt.year == 2020))

sel_vswc_JJA=sel_vswc.sel(time=(uxds_fine2.time.dt.month.isin([6,7,8])&(uxds_fine2.time.dt.year.isin([2020]))))
sel_le_JJA=sel_le.sel(time=(uxds_fine2.time.dt.month.isin([6,7,8])&(uxds_fine2.time.dt.year.isin([2020]))))
sel_h_JJA=sel_h.sel(time=(uxds_fine2.time.dt.month.isin([6,7,8])&(uxds_fine2.time.dt.year.isin([2020]))))
sel_rnet_JJA=sel_rnet.sel(time=(uxds_fine2.time.dt.month.isin([6,7,8])&(uxds_fine2.time.dt.year.isin([2020]))))
sel_t2m_JJA=sel_t2m.sel(time=(uxds_fine2.time.dt.month.isin([6,7,8])&(uxds_fine2.time.dt.year.isin([2020]))))
sel_sp_JJA=sel_sp.sel(time=(uxds_fine2.time.dt.month.isin([6,7,8])&(uxds_fine2.time.dt.year.isin([2020]))))

daily_sel_vswc_DJF=sel_vswc_DJF.resample(time='1D').mean()
daily_sel_le_DJF=sel_le_DJF.resample(time='1D').mean()
daily_sel_h_DJF=sel_h_DJF.resample(time='1D').mean()
daily_sel_rnet_DJF=sel_rnet_DJF.resample(time='1D').mean()
daily_sel_t2m_DJF=sel_t2m_DJF.resample(time='1D').mean()
daily_sel_sp_DJF=sel_sp_DJF.resample(time='1D').mean()

daily_sel_vswc_JJA=sel_vswc_JJA.resample(time='1D').mean()
daily_sel_le_JJA=sel_le_JJA.resample(time='1D').mean()
daily_sel_h_JJA=sel_h_JJA.resample(time='1D').mean()
daily_sel_rnet_JJA=sel_rnet_JJA.resample(time='1D').mean()
daily_sel_t2m_JJA=sel_t2m_JJA.resample(time='1D').mean()
daily_sel_sp_JJA=sel_sp_JJA.resample(time='1D').mean()

daily_sel_le_JJA.shape

(92, 196608)

def saturation_vapor_pressure(T):
    return 6.112*np.exp((17.67*T)/(T+243.5))*0.1  # hPa -> kPa

def delta_svp(T):
    es=saturation_vapor_pressure(T)
    return (4302.645*es)/((T+243.5)**2)  # in kPa/°C

def alpha_piecewise(vswc): # values from literature based on Priestley–Taylor equation
    return xr.where(
        vswc < 0.1, 0.8,
        xr.where(vswc < 0.2, 1.0, 1.26)
    )

def calc_gamma(P_kPa):
    cp=1.013e3      # J/kg/°C
    epsilon=0.622
    lambda_v=2.45e6 # Latent heat of vaporization J/kg
    return (cp*P_kPa)/(epsilon*lambda_v)  # in kPa/°C

def priestley_taylor_ET(vswc, Rnet, T, P_kPa): # assume Rnet=LE+H
    alpha=alpha_piecewise(vswc)
    delta=delta_svp(T)
    gamma=calc_gamma(P_kPa)
    coeff=alpha*(delta/(delta+gamma))
    return coeff*Rnet  # W/m2

def safe_corr(x, y, dim='time', min_count=30):
    valid = x.notnull() & y.notnull()
    count = valid.sum(dim=dim)
    r = xr.corr(x, y, dim=dim)
    return r.where(count >= min_count)
    
def calculate_coupling_index(vswc, Rnet, LE, T, P):
    H=Rnet-LE

    LE_p=priestley_taylor_ET(vswc, Rnet, T, P)
    H_p=Rnet-LE_p
    
    H_xr=H.to_dataset(name='H')
    Hp_xr=H_p.to_dataset(name='Hp')
    T_xr=T.to_dataset(name='T2M')

    H_corr = safe_corr(H, T, dim='time', min_count=30)
    Hp_corr = safe_corr(H_p, T, dim='time', min_count=30)

    PI=H_corr-Hp_corr
    PI.name='coupling_index'

    return PI

%%time
PI_DJF=calculate_coupling_index(daily_sel_vswc_DJF,daily_sel_rnet_DJF,daily_sel_le_DJF,daily_sel_t2m_DJF,daily_sel_sp_DJF).compute()

CPU times: user 2min 23s, sys: 37.5 s, total: 3min
Wall time: 46.3 s

%%time
PI_JJA=calculate_coupling_index(daily_sel_vswc_JJA,daily_sel_rnet_JJA,daily_sel_le_JJA,daily_sel_t2m_JJA,daily_sel_sp_JJA).compute()

CPU times: user 1min 58s, sys: 23.6 s, total: 2min 22s
Wall time: 34 s

PI_JJA.plot()

PI_DJF.plot()

References¶

Miralles, D. G., van den Berg, M. J., Teuling, A. J., & de Jeu, R. A. M. (2012). Soil moisture‐temperature coupling: A multiscale observational analysis. Geophysical Research Letters, 39(21). 10.1029/2012gl053703

Notebooks (work in progress)

Calculate the SM-T coupling using icon_d3hp003

Notebooks (work in progress)

Calculate the LE-prcp correlation