Empty cubes after volume_statistics for certain datasets
**Describe the bug**
CMIP6 EC-Earth ocean 3D data has the whole bottom depth layer set to masked values. Both the variables and the volcello. The volume_statistics preproc computes the weighted mean per layer just fine, and for the bottom layer the results are masked as well. But when it comes to computing the weighted mean over the water volume, the presence of masked values sets the results to NaNs and the data comes out of the preprocessing empty:

`result.append(np.average(column, weights=depth_volume))`

It looks like the function `np.average` is not handling well the presence of layers with only masked values. Would it be correct check and remove this sort of layers before calling `np.average` ?

**Please attach**
I'm working with EC-Earth local data and can't seem to find Ecearth data with volcello on /badc, sorry. 

diff --git a/esmvalcore/preprocessor/_volume.py b/esmvalcore/preprocessor/_volume.py
index d7d5f0203..1e374b6cc 100644
--- a/esmvalcore/preprocessor/_volume.py
+++ b/esmvalcore/preprocessor/_volume.py
@@ -78,13 +78,10 @@ def _create_cube_time(src_cube, data, times):
     Returns
     -------
     cube
-
     .. note::
-
         If there is only one level of interpolation, the resultant cube
         will be collapsed over the associated vertical dimension, and a
         scalar vertical coordinate will be added.
-
     """
     # Get the source cube vertical coordinate and associated dimension.
     src_times = src_cube.coord('time')
@@ -268,12 +265,14 @@ def volume_statistics(
             depth_volume.append(layer_vol)
         # ####
         # Calculate weighted mean over the water volumn
-        result.append(np.average(column, weights=depth_volume))
+        column = np.ma.array(column)
+        depth_volume = np.ma.array(depth_volume)
+        result.append(np.ma.average(column, weights=depth_volume))
 
     # ####
     # Send time series and dummy cube to cube creating tool.
     times = np.array(cube.coord('time').points.astype(float))
-    result = np.array(result)
+    result = np.ma.array(result)
 
     # #####
     # Create a small dummy output array for the output cube
diff --git a/tests/unit/preprocessor/_volume/test_volume.py b/tests/unit/preprocessor/_volume/test_volume.py
index 54b10e188..0e0723160 100644
--- a/tests/unit/preprocessor/_volume/test_volume.py
+++ b/tests/unit/preprocessor/_volume/test_volume.py
@@ -86,7 +86,7 @@ class Test(tests.Test):
     def test_volume_statistics(self):
         """Test to take the volume weighted average of a (2,3,2,2) cube."""
         result = volume_statistics(self.grid_4d, 'mean')
-        expected = np.array([1., 1.])
+        expected = np.ma.array([1., 1.], mask=False)
         self.assert_array_equal(result.data, expected)
 
     def test_volume_statistics_long(self):
@@ -97,7 +97,27 @@ class Test(tests.Test):
         different methods for small and large cubes.
         """
         result = volume_statistics(self.grid_4d_2, 'mean')
-        expected = np.array([1., 1., 1., 1.])
+        expected = np.ma.array([1., 1., 1., 1.], mask=False)
+        self.assert_array_equal(result.data, expected)
+
+    def test_volume_statistics_masked_level(self):
+        """
+        Test to take the volume weighted average of a (2,3,2,2) cube
+        where the last depth level is fully masked.
+        """
+        self.grid_4d.data[:, -1, :, :] = np.ma.masked_all((2, 2, 2))
+        result = volume_statistics(self.grid_4d, 'mean')
+        expected = np.ma.array([1., 1.], mask=False)
+        self.assert_array_equal(result.data, expected)
+
+    def test_volume_statistics_masked_timestep(self):
+        """
+        Test to take the volume weighted average of a (2,3,2,2) cube
+        where the first timestep is fully masked.
+        """
+        self.grid_4d.data[0, :, :, :] = np.ma.masked_all((3, 2, 2))
+        result = volume_statistics(self.grid_4d, 'mean')
+        expected = np.ma.array([1., 1], mask=[True, False])
         self.assert_array_equal(result.data, expected)
 
     def test_depth_integration_1d(self):
['tests/unit/preprocessor/_volume/test_volume.py', 'esmvalcore/preprocessor/_volume.py']