diff --git a/earthdiagnostics/ocean/density.py b/earthdiagnostics/ocean/density.py
index d4d4aea489c0ef4916727764385fc1ce29f0829f..69b070c487f15c6bbe49362ecbcc51d783ca1283 100644
--- a/earthdiagnostics/ocean/density.py
+++ b/earthdiagnostics/ocean/density.py
@@ -134,7 +134,7 @@ class Density(Diagnostic):
                     np.full(bigthetao.shape[1:],
                             ref_pressure, dtype=np.float32)
                 ))
-            sigma_values = np.stack(sigma_values)
+            sigma_values = np.ma.masked_invalid(np.stack(sigma_values))
             sigma_cube = bigthetao.copy(sigma_values)
             sigma_cube.var_name = f'sigma{sigma}'
             sigma_cube.standard_name = 'sea_water_sigma_theta'
@@ -144,7 +144,7 @@ class Density(Diagnostic):
             )
             sigma_cube.units = 'kg m-3'
             temp = TempFile.get()
-            iris.save(sigma_cube, temp, zlib=True)
+            iris.save(sigma_cube, temp, zlib=True, fill_value=1e20)
             del sigma_cube
             del sigma_values
             self.sigma[sigma].set_local_file(temp)