diff --git a/lib/iris/_lazy_data.py b/lib/iris/_lazy_data.py
index 985e79af13..ba16dc4a8a 100644
--- a/lib/iris/_lazy_data.py
+++ b/lib/iris/_lazy_data.py
@@ -126,6 +126,34 @@ def as_concrete_data(data):
     return data
 
 
+def lazy_masked_fill_value(data):
+    """
+    Return the fill value of a lazy masked array.
+
+    Args:
+
+    * data:
+        A dask array, NumPy `ndarray` or masked array
+
+    Returns:
+        The fill value of `data` if `data` represents a masked array, or None.
+
+    """
+    # If lazy, get the smallest slice of the data from which we can retrieve
+    # the fill_value.
+    if is_lazy_data(data):
+        inds = tuple([0] * (data.ndim-1))
+        smallest_slice = data[inds][:0]
+        data = as_concrete_data(smallest_slice)
+
+    # Now get the fill value.
+    fill_value = None
+    if ma.isMaskedArray(data):
+        fill_value = data.fill_value
+
+    return fill_value
+
+
 def multidim_lazy_stack(stack):
     """
     Recursively build a multidimensional stacked dask array.
diff --git a/lib/iris/fileformats/netcdf.py b/lib/iris/fileformats/netcdf.py
index e411d3ac1b..a58a583c67 100644
--- a/lib/iris/fileformats/netcdf.py
+++ b/lib/iris/fileformats/netcdf.py
@@ -56,7 +56,7 @@
 import iris.fileformats._pyke_rules
 import iris.io
 import iris.util
-from iris._lazy_data import as_lazy_data
+from iris._lazy_data import as_lazy_data, lazy_masked_fill_value
 
 # Show Pyke inference engine statistics.
 DEBUG = False
@@ -1946,12 +1946,16 @@ def set_packing_ncattrs(cfvar):
 
         else:
             # Create the cube CF-netCDF data variable.
-            # Explicitly assign the fill_value, which will be the type default
-            # in the case of an unmasked array.
+            # Set `fill_value` if the data array is masked. If the data array
+            # is lazy masked, we realise the smallest possible slice of the
+            # array and retrieve the fill value from that.
             if packing is None:
-                fill_value = None
                 if not cube.has_lazy_data() and ma.isMaskedArray(cube.data):
                     fill_value = cube.data.fill_value
+                elif cube.has_lazy_data():
+                    fill_value = lazy_masked_fill_value(cube.lazy_data())
+                else:
+                    fill_value = None
                 dtype = cube.dtype.newbyteorder('=')
 
             cf_var = self._dataset.createVariable(
diff --git a/lib/iris/tests/unit/lazy_data/test_lazy_masked_fill_value.py b/lib/iris/tests/unit/lazy_data/test_lazy_masked_fill_value.py
new file mode 100644
index 0000000000..aaaefea657
--- /dev/null
+++ b/lib/iris/tests/unit/lazy_data/test_lazy_masked_fill_value.py
@@ -0,0 +1,86 @@
+# (C) British Crown Copyright 2017, Met Office
+#
+# This file is part of Iris.
+#
+# Iris is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the
+# Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# Iris is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with Iris.  If not, see <http://www.gnu.org/licenses/>.
+"""Test the function :func:`iris._lazy data.lazy_masked_fill_value`."""
+
+from __future__ import (absolute_import, division, print_function)
+from six.moves import (filter, input, map, range, zip)  # noqa
+
+# Import iris.tests first so that some things can be initialised before
+# importing anything else.
+import iris.tests as tests
+
+import dask.array as da
+import numpy as np
+import numpy.ma as ma
+
+from iris._lazy_data import lazy_masked_fill_value, _MAX_CHUNK_SIZE
+
+
+class Test_as_lazy_data(tests.IrisTest):
+    def setUp(self):
+        shape = (2, 3, 4)
+        data = np.arange(24).reshape(shape)
+        mask = np.zeros(shape)
+        mask[data % 5 == 1] = 1
+        self.fill_value = 9999
+        self.m = ma.masked_array(data, mask=mask, fill_value=self.fill_value)
+        self.dm = da.from_array(self.m, asarray=False,
+                                chunks=_MAX_CHUNK_SIZE)
+
+    def test_lazy_masked_ND(self):
+        fill_value = lazy_masked_fill_value(self.dm)
+        self.assertEqual(fill_value, self.fill_value)
+
+    def test_lazy_masked_0D(self):
+        data = self.dm[0, 0, :1]
+        fill_value = lazy_masked_fill_value(data)
+        self.assertEqual(fill_value, self.fill_value)
+
+    def test_lazy_masked_1D(self):
+        data = self.dm[0, 0, :]
+        fill_value = lazy_masked_fill_value(data)
+        self.assertEqual(fill_value, self.fill_value)
+
+    def test_lazy_masked_2D(self):
+        data = self.dm[0, :]
+        fill_value = lazy_masked_fill_value(data)
+        self.assertEqual(fill_value, self.fill_value)
+
+    def test_real_masked(self):
+        fill_value = lazy_masked_fill_value(self.m)
+        self.assertEqual(fill_value, self.fill_value)
+
+    def test_lazy_unmasked(self):
+        data = da.from_array(self.m.filled(),
+                             chunks=_MAX_CHUNK_SIZE)
+        fill_value = lazy_masked_fill_value(data)
+        self.assertIsNone(fill_value)
+
+    def test_real_unmasked(self):
+        data = self.m.filled()
+        fill_value = lazy_masked_fill_value(data)
+        self.assertIsNone(fill_value)
+
+    def test_data_not_realised(self):
+        # Check that only the zero-element slice is realised.
+        data = self.dm[0, :]
+        lazy_masked_fill_value(data)
+        self.assertIsInstance(data, da.core.Array)
+
+
+if __name__ == '__main__':
+    tests.main()