WIP: Introduce a second version of the Holland 2010 model

climada/hazard/test/test_trop_cyclone.py

@@ -33,8 +33,8 @@
 from climada.hazard.tc_tracks import TCTracks
 from climada.hazard.trop_cyclone import (
     TropCyclone, get_close_centroids, _vtrans, _B_holland_1980, _bs_holland_2008,
-    _v_max_s_holland_2008, _x_holland_2010, _stat_holland_1980, _stat_holland_2010,
-    _stat_er_2011, tctrack_to_si, MBAR_TO_PA, KM_TO_M, H_TO_S,
+    _bs_holland_2010_v2, _v_max_s_holland_2008, _x_holland_2010, _stat_holland_1980,
+    _stat_holland_2010, _stat_er_2011, tctrack_to_si, MBAR_TO_PA, KM_TO_M, H_TO_S,
 )
 from climada.hazard.centroids.centr import Centroids
 import climada.hazard.test as hazard_test
@@ -139,19 +139,23 @@ def test_windfield_models(self):
                 24.745521, 25.596484, 26.475329, 24.690914, 28.650107, 31.584395,
                 21.723546, 26.140293, 28.94964,  28.051915, 18.49378, 35.312152,
             ],
-            # Holland 1980 and Emanuel & Rotunno 2011 use recorded wind speeds, while the above use
-            # pressure values only. That's why the results are so different:
+            # Holland 1980, version 2 of Holland 2010 and Emanuel & Rotunno 2011 use recorded wind speeds,
+            # while the above use pressure values only. That's why the results are so different:
             "H1980": [21.376807, 21.957217, 22.569568, 21.284351, 24.254226, 26.971303,
                       19.220149, 21.984516, 24.196388, 23.449116,  0, 31.550207],
             "ER11": [23.565332, 24.931413, 26.360758, 23.490333, 29.601171, 34.522795,
                      18.996389, 26.102109, 30.780737, 29.498453,  0, 38.368805],
+            "H10_v2": [
+                28.067253, 28.544574, 28.975862, 27.907048, 30.450187, 32.091651,
+                25.79227 , 28.140679, 29.615011, 28.693995, 22.191709, 32.245298
+            ],
         }
 
         tc_track = TCTracks.from_processed_ibtracs_csv(TEST_TRACK)
         tc_track.equal_timestep()
         tc_track.data = tc_track.data[:1]
 
-        for model in ["H08", "H10", "H1980", "ER11"]:
+        for model in ["H08", "H10", "H1980", "ER11", "H10_v2"]:
             tc_haz = TropCyclone.from_tracks(tc_track, centroids=CENTR_TEST_BRB, model=model)
             np.testing.assert_array_almost_equal(
                 tc_haz.intensity[0, intensity_idx].toarray()[0], intensity_values[model])
@@ -299,6 +303,17 @@ def test_bs_holland_2008_pass(self):
         np.testing.assert_array_almost_equal(
             si_track["hol_b"], [np.nan, 1.27085617, 1.26555341])
 
+    def test_bs_holland_2010_v2_pass(self):
+        """Test _bs_holland_2010_v2 function. Compare to MATLAB reference."""
+        si_track = xr.Dataset({
+            "env": ("time", MBAR_TO_PA * np.array([1010, 1010, 1010, 1010])),
+            "cen": ("time", MBAR_TO_PA * np.array([1005.2585, 1005.2633, 1005.2682, 1005.2731])),
+            "vmax": ("time",  [np.nan, 22.5, 25.4, 42.5]),
+        })
+        _bs_holland_2010_v2(si_track)
+        np.testing.assert_array_almost_equal(
+            si_track["hol_b"], [np.nan, 1.75439, 2.238092, 2.5])
+
     def test_v_max_s_holland_2008_pass(self):
         """Test _v_max_s_holland_2008 function."""
         # Numbers analogous to test_B_holland_1980_pass

climada/hazard/trop_cyclone.py

@@ -59,11 +59,13 @@
 DEF_MAX_MEMORY_GB = 8
 """Default value of the memory limit (in GB) for windfield computations (in each thread)."""
 
-MODEL_VANG = {'H08': 0, 'H1980': 1, 'H10': 2, 'ER11': 3}
+MODEL_VANG = {'H08': 0, 'H1980': 1, 'H10': 2, 'ER11': 3, 'H10_v2': 4}
 """Enumerate different symmetric wind field models."""
 
 RHO_AIR = 1.15
 """Air density. Assumed constant, following Holland 1980."""
+RHO_AIR_SURFACE = 1.2
+"""Air density at surface level. Assumed constant, following Holland 2010"""
 
 GRADIENT_LEVEL_TO_SURFACE_WINDS = 0.9
 """Gradient-to-surface wind reduction factor according to the 90%-rule:
@@ -969,7 +971,7 @@
     v_trans_corr[close_centr_msk] = np.fmin(
         1, t_rad_bc[close_centr_msk] / d_centr[close_centr_msk])
 
-    if model in [MODEL_VANG['H08'], MODEL_VANG['H10']]:
+    if model in [MODEL_VANG['H08'], MODEL_VANG['H10'], MODEL_VANG['H10_v2']]:
         # In these models, v_ang_norm already contains vtrans_norm, so subtract it first, before
         # converting to vectors and then adding (vectorial) vtrans again. Make sure to apply the
         # "absorbing factor" in both steps:
@@ -1106,7 +1108,8 @@
         _B_holland_1980(si_track)
     elif model in [MODEL_VANG['H08'], MODEL_VANG['H10']]:
         _bs_holland_2008(si_track)
-
+    elif model in [MODEL_VANG['H10_v2']]:
+        _bs_holland_2010_v2(si_track)
     if model in [MODEL_VANG['H1980'], MODEL_VANG['H08']]:
         result = _stat_holland_1980(
             si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic,
@@ -1120,6 +1123,9 @@
         result = _stat_holland_2010(si_track, d_centr, close_centr_msk, hol_x)
     elif model == MODEL_VANG['ER11']:
         result = _stat_er_2011(si_track, d_centr, close_centr_msk, cyclostrophic=cyclostrophic)
+    elif model in [MODEL_VANG['H10_v2'],]:
+        hol_x = _x_holland_2010(si_track, d_centr, close_centr_msk)
+        result = _stat_holland_2010(si_track, d_centr, close_centr_msk, hol_x)
     else:
         raise NotImplementedError
 
@@ -1255,6 +1261,8 @@
 
 def _bs_holland_2008(si_track: xr.Dataset):
     """Holland's 2008 b-value estimate for sustained surface winds.
+    (This is also one option of how to estimate the b-value in Holland 2010,
+    for the other option consult '_bs_holland_2010_v2'.)
 
     The result is stored in place as a new data variable "hol_b".
 
@@ -1311,6 +1319,44 @@
     )
     si_track["hol_b"] = np.clip(si_track["hol_b"], 1, 2.5)
 
+def _bs_holland_2010_v2(si_track: xr.Dataset):
+    """Holland 2010's second version of how to estimate b-value  for sustained surface winds.
+    For version 1 see "_bs_holland_2008"
+
+    The result is stored in place as a new data variable "hol_b".
+
+    Like the original 1980 formula (see `_B_holland_1980`), this approach does also require
+    wind speed measurements, if the wind speed measurements are not available or not  reliable, consider the second
+    option proposed in Holland 2010 (and Holland 2008) implemented as "_bs_holland_2008"
+
+    The parameter applies to 1-minute sustained winds at 10 meters above ground.
+    It is taken from equation (7) in the following paper:
+
+    Holland et al. (2010): A Revised Model for Radial Profiles of Hurricane Winds. Monthly
+    Weather Review 138(12): 4393–4401. https://doi.org/10.1175/2010MWR3317.1
+
+    For reference, it reads
+    b_s = vmax^2 * rho^e / ( 100 * (penv - pcen) )
+    where:
+        rho is the air density ρ (in kg m−3) at the surface level
+        !penv and pcen is assumed to be in hPa in this formula - not Pa, as in our tracks
+
+    Parameters
+    ----------
+    si_track : xr.Dataset
+        Output of `tctrack_to_si`. The data variables used by this function are
+        "cen",  "env", and "vmax". The result is stored in place as a new data
+        variable "hol_b".
+    """
+
+    si_track["hol_b"] = (
+            si_track["vmax"]**2 * RHO_AIR_SURFACE**np.exp(1)
+            / ( si_track["env"] - si_track["cen"] ) # we do not need the factor 100 as in the original
+        # formula, because environmental pressure and central pressure are already saved in Pa, not hPa
+    )
+    si_track["hol_b"] = np.clip(si_track["hol_b"], 0.4, 2.5)  # reconsider the clip as b_s is b_g/1.6, but does this also relate to limits?
+
+
 def _v_max_s_holland_2008(si_track: xr.Dataset):
     """Compute maximum surface winds from pressure according to Holland 2008.
 
@@ -1433,7 +1479,7 @@
 
     # compute peripheral exponent from second measurement
     r_max_norm = (r_max / r_n)**hol_b
-    x_n = np.log(v_n / v_max_s) / np.log(r_max_norm * np.exp(1 - r_max_norm))
+    x_n = np.log(v_n / v_max_s) / np.log(r_max_norm * np.exp(1 - r_max_norm))  # holland 2010 equation 6 solved for x
 
     # linearly interpolate between max exponent and peripheral exponent
     x_max = 0.5