diff --git a/reproject/mosaicking/coadd.py b/reproject/mosaicking/coadd.py
index e9c49d3c1..7a99e2a9c 100644
--- a/reproject/mosaicking/coadd.py
+++ b/reproject/mosaicking/coadd.py
@@ -276,6 +276,7 @@ def reproject_and_coadd(
         if combine_function == "mean":
             with np.errstate(invalid="ignore"):
                 output_array /= output_footprint
+                output_array[output_footprint == 0] = 0
 
     elif combine_function in ("first", "last", "min", "max"):
         for array in arrays:
diff --git a/reproject/mosaicking/tests/reference/test_coadd_solar_map.fits b/reproject/mosaicking/tests/reference/test_coadd_solar_map.fits
index 1ec779a00..7d8747d26 100644
Binary files a/reproject/mosaicking/tests/reference/test_coadd_solar_map.fits and b/reproject/mosaicking/tests/reference/test_coadd_solar_map.fits differ
diff --git a/reproject/mosaicking/tests/test_coadd.py b/reproject/mosaicking/tests/test_coadd.py
index fb3f757a8..f259e12fd 100644
--- a/reproject/mosaicking/tests/test_coadd.py
+++ b/reproject/mosaicking/tests/test_coadd.py
@@ -210,7 +210,7 @@ def test_coadd_background_matching(self, reproject_function):
 
         assert_allclose(array - np.mean(array), self.array - np.mean(self.array), atol=ATOL)
 
-    @pytest.mark.parametrize("combine_function", ["first", "last", "min", "max", "sum"])
+    @pytest.mark.parametrize("combine_function", ["first", "last", "min", "max", "sum", "mean"])
     @pytest.mark.parametrize("match_background", [True, False])
     def test_footprint_correct(self, reproject_function, combine_function, match_background):
         # Test that the output array is zero outside the returned footprint