landfast (or static) sea ice calculation

[dpath2o]

I would like to provide a function to cosima-recipes that allows one to compute landfast sea ice.

Landfast sea ice is defined as sea ice that is fixed to the coast, shelf, tongue, iceberg and does not move within a time period. In a numerical model that reports sea ice at daily time intervals then one needs to define when the sea ice has become fastened. This is not as straightforward as finding when sea ice velocities are zero and this due to primarily wind shear (and likely in the future surface gravity waves) that are imparting dynamic forces that are translating to even very small velocities. Hence a two-week window is used to define landfast ice as a shorter period could lead to false assessments for low wind events and a longer period would provide too much temporal smoothing. Lemieux 2015 provides further (a more detailed) justification for averaging window, with another important point that the 'fast ice threshold' (5e-4 m/s) needs to be larger than the basal stress parameter (CICE default 5e-e5 m/s).

With the above being said, it is noted that ACCESS-OM2 does not include fast ice dynamics in the sea ice model. Hence when computing/reporting/masking for fast ice in ACCESS-OM2, one is cautioned in calling it 'fast ice' and potentially a more appropriate name for such simulated sea ice could be 'static ice'.

[dpath2o]

Here are the functions that I'll be adding.

def mask_hemisphere(DS, var_name=None, grid_type='t', hemisphere='south'):
        """
        Apply a hemisphere mask to a dataset variable; *not* required for NSIDC data

        Parameters:
        -----------
        DS : xarray.Dataset
            The dataset containing the variable to be masked.
        var_name : str
            The name of the variable within the dataset to be masked.
        grid_type : str, optional
            The type of grid to use ('t', 'u', 'ot', 'ou'). Default is 't'.
        hemisphere : str, optional
            The hemisphere to mask ('north' or 'south'). Default is 'south'.

        Returns:
        --------
        xarray.DataArray
            The masked DataArray for the specified variable.

        Raises:
        -------
        ValueError
            If the hemisphere is not 'north' or 'south'.

        Example:
        --------
        masked_da = analysis.mask_hemisphere(DS, var_name='aice', grid_type='t', hemisphere='north')
        """
        DA = DS[var_name]
        if hemisphere in ['southern', 'sh', 'south', 's', 'sthn', 'sth', 'SH']:
            hemisphere = 'south'
        elif hemisphere in ['northern', 'nh', 'north', 'n', 'nthn', 'nth', 'NH']:
            hemisphere = 'north'
        else:
            raise ValueError("Hemisphere must be 'north' or 'south'.")
        op = (lambda x: x < 0) if hemisphere == 'south' else (lambda x: x > 0)
        mask = op(DS[self.grid_shorthand[grid_type]])
        return DA.where(mask)

def compute_fast_ice(ds, hemisphere='SH', fic=False, fia=False, fie=False, aice_thresh=0.15, fi_thesh=5e-4):
    """
    Compute landfast sea ice concentration and/or area and/or extent.

    Parameters:
    ds (xarray.Dataset): The dataset containing the necessary variables.
    hemisphere (str): 'SH' for Southern Hemisphere or 'NH' for Northern Hemisphere. Default is 'SH'.
     fic (bool): If True, compute and return static ice concentration.
     fia (bool): If True, compute and return static ice area.
     fie (bool): If True, compute and return static ice extent.
     aice_thresh (float): a value that is used for determining if a grid cell has a sufficient concentration of sea ice to be defined as a grid cell containing sea ice.
     fi_thresh (float): a value that defines the upper limit of the ice speed in m/s, everthing less than this value will be considered fast ice. 

    Returns:
    tuple: (xarray.DataArray, xarray.DataArray, xarray.DataArray)
           fast ice concentration, fast ice area, and fast ice extent depending on the parameters.
    """
    # Helper function to mask and resample variables
    def mask_and_resample(var_name, grid_type, hemisphere, ds):
        return mask_hemisphere(ds, var_name=var_name, grid_type=grid_type, hemisphere=hemisphere).sortby('time').resample(time='2W').mean()
    # Load and process variables
    uvel = mask_and_resample('uvel', 'u', hemisphere, ds)
    vvel = mask_and_resample('vvel', 'u', hemisphere, ds)
    aice = mask_and_resample('aice', 't', hemisphere, ds)
    # Compute ice speed
    ispd = np.sqrt(uvel**2 + vvel**2)
    # Create fast ice mask
    mask_stic = (aice > aice_thresh) & (ispd < fi_thresh)
    fic_result = None
    fia_result = None
    fie_result = None
    # Compute static ice concentration    
    if fic:
        fic_result = aice.where(mask_stic)
    # Compute ice area
    if fia:
        area = mask_hemisphere(ds, var_name='uarea', grid_type='u', hemisphere=hemisphere)
        fia_result = area.where(mask_stic).sum(dim=['nj', 'ni']) / 1e9 #1e3 km^2
    # Compute ice extent
    if fie:
        area = mask_hemisphere(ds, var_name='uarea', grid_type='u', hemisphere=hemisphere)
        fie_result = mask_stic.sum(dim=['nj', 'ni'], skipna=True) / 1e9  #1e3 km^2
    # return results
    return fic_result, fia_result, fie_result

[anton-seaice]

Thanks Dan

Mostly cosima-recipes uses a minimal number of functions. It would also be neat to use cf compliance.

i.e. the SH mask could be

data_array.cf.sel(latitude=slice(-90,0))

and

.sum(dim=['nj','ni']) could be .cf.sum(dim='longitude','latitude')