diff --git a/src/xclim/indices/_threshold.py b/src/xclim/indices/_threshold.py
index 4404461c7..7fac2949c 100644
--- a/src/xclim/indices/_threshold.py
+++ b/src/xclim/indices/_threshold.py
@@ -22,7 +22,9 @@
 )
 from xclim.indices import run_length as rl
 from xclim.indices.generic import (
+    bivariate_count_occurrences,
     compare,
+    count_occurrences,
     cumulative_difference,
     domain_count,
     first_day_threshold_reached,
@@ -3457,7 +3459,7 @@ def wet_spell_frequency(
         Resampling frequency.
     resample_before_rl : bool
         Determines if the resampling should take place before or after the run length encoding (or a similar algorithm) is applied to runs.
-    op : {"sum","min", "max", "mean"}
+    op : {"sum", "min", "max", "mean"}
         Operation to perform on the window.
         Default is "sum", which checks that the sum of accumulated precipitation over the whole window is more than the
         threshold.
@@ -3647,14 +3649,15 @@ def wet_spell_max_length(
 def holiday_snow_days(
     snd: xarray.DataArray,
     snd_thresh: Quantified = "20 mm",
+    op: str = ">=",
     date_start: str = "12-25",
     date_end: str | None = None,
     freq: str = "YS",
 ) -> xarray.DataArray:  # numpydoc ignore=SS05
-    r"""
+    """
     Christmas Days.
 
-    Whether there is a significant amount of snow on the ground on December 25th (or around that time).
+    Whether there is a significant amount of snow on the ground on December 25th (or a given date range).
 
     Parameters
     ----------
@@ -3662,6 +3665,8 @@ def holiday_snow_days(
         Surface snow depth.
     snd_thresh : Quantified
         Threshold snow amount. Default: 20 mm.
+    op : {">", "gt", ">=", "ge"}
+        Comparison operation. Default: ">=".
     date_start : str
         Beginning of analysis period. Default: "12-25" (December 25th).
     date_end : str, optional
@@ -3680,15 +3685,15 @@ def holiday_snow_days(
     ----------
     https://www.canada.ca/en/environment-climate-change/services/weather-general-tools-resources/historical-christmas-snowfall-data.html
     """
-    snow_depth = convert_units_to(snd, "m")
-    snow_depth_thresh = convert_units_to(snd_thresh, "m")
-    snow_depth_constrained = select_time(
-        snow_depth,
+    snd_constrained = select_time(
+        snd,
         drop=True,
         date_bounds=(date_start, date_start if date_end is None else date_end),
     )
 
-    xmas_days = (snow_depth_constrained >= snow_depth_thresh).resample(time=freq).sum()
+    xmas_days = count_occurrences(
+        snd_constrained, snd_thresh, freq, op, constrain=[">=", ">"]
+    )
 
     xmas_days = xmas_days.assign_attrs({"units": "days"})
     return xmas_days
@@ -3705,6 +3710,8 @@ def holiday_snow_and_snowfall_days(
     prsn: xarray.DataArray | None = None,
     snd_thresh: Quantified = "20 mm",
     prsn_thresh: Quantified = "1 cm",
+    snd_op: str = ">=",
+    prsn_op: str = ">=",
     date_start: str = "12-25",
     date_end: str | None = None,
     freq: str = "YS-JUL",
@@ -3724,6 +3731,10 @@ def holiday_snow_and_snowfall_days(
         Threshold snow amount. Default: 20 mm.
     prsn_thresh : Quantified
         Threshold snowfall flux. Default: 1 mm.
+    snd_op : {">", "gt", ">=", "ge"}
+        Comparison operation for snow depth. Default: ">=".
+    prsn_op : {">", "gt", ">=", "ge"}
+        Comparison operation for snowfall flux. Default: ">=".
     date_start : str
         Beginning of analysis period. Default: "12-25" (December 25th).
     date_end : str, optional
@@ -3742,31 +3753,32 @@ def holiday_snow_and_snowfall_days(
     ----------
     https://www.canada.ca/en/environment-climate-change/services/weather-general-tools-resources/historical-christmas-snowfall-data.html
     """
-    snow_depth = convert_units_to(snd, "m")
-    snow_depth_thresh = convert_units_to(snd_thresh, "m")
-    snow_depth_constrained = select_time(
-        snow_depth,
+    snd_constrained = select_time(
+        snd,
         drop=True,
         date_bounds=(date_start, date_start if date_end is None else date_end),
     )
 
-    snowfall_rate = rate2amount(
+    prsn_mm = rate2amount(
         convert_units_to(prsn, "mm day-1", context="hydro"), out_units="mm"
     )
-    snowfall_thresh = convert_units_to(prsn_thresh, "mm", context="hydro")
-    snowfall_constrained = select_time(
-        snowfall_rate,
+    prsn_mm_constrained = select_time(
+        prsn_mm,
         drop=True,
         date_bounds=(date_start, date_start if date_end is None else date_end),
     )
 
-    perfect_xmas_days = (
-        (
-            (snow_depth_constrained >= snow_depth_thresh)
-            & (snowfall_constrained >= snowfall_thresh)
-        )
-        .resample(time=freq)
-        .sum()
+    perfect_xmas_days = bivariate_count_occurrences(
+        data_var1=snd_constrained,
+        data_var2=prsn_mm_constrained,
+        threshold_var1=snd_thresh,
+        threshold_var2=prsn_thresh,
+        op_var1=snd_op,
+        op_var2=prsn_op,
+        freq=freq,
+        var_reducer="all",
+        constrain_var1=[">=", ">"],
+        constrain_var2=[">=", ">"],
     )
 
     perfect_xmas_days = perfect_xmas_days.assign_attrs({"units": "days"})
diff --git a/src/xclim/indices/generic.py b/src/xclim/indices/generic.py
index e9dd96920..f8ead5977 100644
--- a/src/xclim/indices/generic.py
+++ b/src/xclim/indices/generic.py
@@ -34,6 +34,7 @@
 __all__ = [
     "aggregate_between_dates",
     "binary_ops",
+    "bivariate_count_occurrences",
     "bivariate_spell_length_statistics",
     "compare",
     "count_level_crossings",
@@ -903,9 +904,9 @@ def count_occurrences(
     """
     Calculate the number of times some condition is met.
 
-    First, the threshold is transformed to the same standard_name and units as the input data:
+    First, the threshold is transformed to the same standard_name and units as the input data;
     Then the thresholding is performed as condition(data, threshold),
-    i.e. if condition is `<`, then this counts the number of times `data < threshold`:
+    i.e. if condition is `<`, then this counts the number of times `data < threshold`;
     Finally, count the number of occurrences when condition is met.
 
     Parameters
@@ -934,6 +935,79 @@ def count_occurrences(
     return to_agg_units(out, data, "count", dim="time")
 
 
+@declare_relative_units(threshold_var1="<data_var1>", threshold_var2="<data_var2>")
+def bivariate_count_occurrences(
+    *,
+    data_var1: xr.DataArray,
+    data_var2: xr.DataArray,
+    threshold_var1: Quantified,
+    threshold_var2: Quantified,
+    freq: str,
+    op_var1: str,
+    op_var2: str,
+    var_reducer: str,
+    constrain_var1: Sequence[str] | None = None,
+    constrain_var2: Sequence[str] | None = None,
+) -> xr.DataArray:
+    """
+    Calculate the number of times some conditions are met for two variables.
+
+    First, the thresholds are transformed to the same standard_name and units as their corresponding input data;
+    Then the thresholding is performed as condition(data, threshold) for each variable,
+    i.e. if condition is `<`, then this counts the number of times `data < threshold`;
+    Then the conditions are combined according to `var_reducer`;
+    Finally, the number of occurrences where conditions are met for "all" or "any" events are counted.
+
+    Parameters
+    ----------
+    data_var1 : xr.DataArray
+        An array.
+    data_var2 : xr.DataArray
+        An array.
+    threshold_var1 : Quantified
+        Threshold for data variable 1.
+    threshold_var2 : Quantified
+        Threshold for data variable 2.
+    freq : str
+        Resampling frequency defining the periods as defined in :ref:`timeseries.resampling`.
+    op_var1 : {">", "gt", "<", "lt", ">=", "ge", "<=", "le", "==", "eq", "!=", "ne"}
+        Logical operator for data variable 1. e.g. arr > thresh.
+    op_var2 : {">", "gt", "<", "lt", ">=", "ge", "<=", "le", "==", "eq", "!=", "ne"}
+        Logical operator for data variable 2. e.g. arr > thresh.
+    var_reducer : {"all", "any"}
+        The condition must either be fulfilled on *all* or *any* variables for the period to be considered an occurrence.
+    constrain_var1 : sequence of str, optional
+        Optionally allowed comparison operators for variable 1.
+    constrain_var2 : sequence of str, optional
+        Optionally allowed comparison operators for variable 2.
+
+    Returns
+    -------
+    xr.DataArray
+        The DataArray of counted occurrences.
+
+    Notes
+    -----
+    Sampling and variable units are derived from `data_var1`.
+    """
+    threshold_var1 = convert_units_to(threshold_var1, data_var1)
+    threshold_var2 = convert_units_to(threshold_var2, data_var2)
+
+    cond_var1 = compare(data_var1, op_var1, threshold_var1, constrain_var1)
+    cond_var2 = compare(data_var2, op_var2, threshold_var2, constrain_var2)
+
+    if var_reducer == "all":
+        cond = cond_var1 & cond_var2
+    elif var_reducer == "any":
+        cond = cond_var1 | cond_var2
+    else:
+        raise ValueError(f"Unsupported value for var_reducer: {var_reducer}")
+
+    out = cond.resample(time=freq).sum()
+
+    return to_agg_units(out, data_var1, "count", dim="time")
+
+
 def diurnal_temperature_range(
     low_data: xr.DataArray, high_data: xr.DataArray, reducer: str, freq: str
 ) -> xr.DataArray: