diff --git a/modules/util/util.py b/modules/util/util.py
index e8064c4d13ed3700240e9b83c145c870bd2b2633..e049667eed646fd88cb27e728cd86b94943fdd83 100644
--- a/modules/util/util.py
+++ b/modules/util/util.py
@@ -608,9 +608,13 @@ def write_icing_file(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y, lon
     dim_0_name = 'pixel_elements'
     dim_1_name = 'scan_lines'
 
+    prob_s = []
+    pred_s = []
+
     flt_lvls = list(preds_dct.keys())
     for flvl in flt_lvls:
         preds = preds_dct[flvl]
+        pred_s.append(preds)
         icing_pred_ds = h5f_out.create_dataset('icing_prediction_level_'+flt_level_ranges_str[flvl], data=preds, dtype='i2')
         icing_pred_ds.attrs.create('coordinates', data='y x')
         icing_pred_ds.attrs.create('grid_mapping', data='Projection')
@@ -620,6 +624,7 @@ def write_icing_file(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y, lon
 
     for flvl in flt_lvls:
         probs = probs_dct[flvl]
+        prob_s.append(probs)
         icing_prob_ds = h5f_out.create_dataset('icing_probability_level_'+flt_level_ranges_str[flvl], data=probs, dtype='f4')
         icing_prob_ds.attrs.create('coordinates', data='y x')
         icing_prob_ds.attrs.create('grid_mapping', data='Projection')
@@ -627,6 +632,27 @@ def write_icing_file(clvrx_str_time, output_dir, preds_dct, probs_dct, x, y, lon
         icing_prob_ds.dims[0].label = dim_0_name
         icing_prob_ds.dims[1].label = dim_1_name
 
+    prob_s = np.stack(prob_s, axis=-1)
+    max_prob = np.max(prob_s, axis=2)
+
+    icing_prob_ds = h5f_out.create_dataset('max_icing_probability_column', data=max_prob, dtype='f4')
+    icing_prob_ds.attrs.create('coordinates', data='y x')
+    icing_prob_ds.attrs.create('grid_mapping', data='Projection')
+    icing_prob_ds.attrs.create('missing', data=-1.0)
+    icing_prob_ds.dims[0].label = dim_0_name
+    icing_prob_ds.dims[1].label = dim_1_name
+
+    pred_s = np.stack(pred_s, axis=-1)
+    sum_pred = np.sum(pred_s, axis=2)
+    sum_pred = np.where(sum_pred < 0, -1, sum_pred)
+
+    icing_pred_ds = h5f_out.create_dataset('sum_icing_column', data=sum_pred, dtype='i2')
+    icing_pred_ds.attrs.create('coordinates', data='y x')
+    icing_pred_ds.attrs.create('grid_mapping', data='Projection')
+    icing_pred_ds.attrs.create('missing', data=-1)
+    icing_pred_ds.dims[0].label = dim_0_name
+    icing_pred_ds.dims[1].label = dim_1_name
+
     lon_ds = h5f_out.create_dataset('longitude', data=lons, dtype='f4')
     lon_ds.attrs.create('units', data='degrees_east')
     lon_ds.attrs.create('long_name', data='icing prediction longitude')