diff --git a/docs/src/whatsnew/latest.rst b/docs/src/whatsnew/latest.rst
index 5c62a27830..e4906894ee 100644
--- a/docs/src/whatsnew/latest.rst
+++ b/docs/src/whatsnew/latest.rst
@@ -128,6 +128,10 @@ This document explains the changes made to Iris for this release
    code comments, this was supposedly already the case, but there were several bugs
    and loopholes. (:issue:`1897`, :pull:`4767`)
 
+#. `@rcomer`_ modified cube arithmetic to handle mismatches in the cube's data
+   array type.  This prevents masks being lost in some cases and therefore
+   resolves :issue:`2987`.  (:pull:`3790`)
+
 
 💣 Incompatible Changes
 =======================
diff --git a/lib/iris/analysis/maths.py b/lib/iris/analysis/maths.py
index 1cbc90cc60..94f01c5140 100644
--- a/lib/iris/analysis/maths.py
+++ b/lib/iris/analysis/maths.py
@@ -837,6 +837,20 @@ def unary_func(lhs):
             raise TypeError(emsg)
         return data
 
+    if in_place and not cube.has_lazy_data():
+        # In-place arithmetic doesn't work if array type of LHS is less complex
+        # than RHS.
+        if iris._lazy_data.is_lazy_data(rhs):
+            cube.data = cube.lazy_data()
+        elif ma.is_masked(rhs) and not isinstance(cube.data, ma.MaskedArray):
+            cube.data = ma.array(cube.data)
+
+    elif isinstance(
+        cube.core_data(), ma.MaskedArray
+    ) and iris._lazy_data.is_lazy_data(rhs):
+        # Workaround for #2987.  numpy#15200 discusses the general problem.
+        cube = cube.copy(cube.lazy_data())
+
     result = _math_op_common(
         cube,
         unary_func,
diff --git a/lib/iris/tests/test_basic_maths.py b/lib/iris/tests/test_basic_maths.py
index 24f2b89442..6c08dc1f9e 100644
--- a/lib/iris/tests/test_basic_maths.py
+++ b/lib/iris/tests/test_basic_maths.py
@@ -687,12 +687,12 @@ def setUp(self):
         self.data_1u = np.array([[9, 9, 9], [8, 8, 8]], dtype=np.uint64)
         self.data_2u = np.array([[3, 3, 3], [2, 2, 2]], dtype=np.uint64)
 
-        self.cube_1f = Cube(self.data_1f)
-        self.cube_2f = Cube(self.data_2f)
-        self.cube_1i = Cube(self.data_1i)
-        self.cube_2i = Cube(self.data_2i)
-        self.cube_1u = Cube(self.data_1u)
-        self.cube_2u = Cube(self.data_2u)
+        self.cube_1f = Cube(self.data_1f.copy())
+        self.cube_2f = Cube(self.data_2f.copy())
+        self.cube_1i = Cube(self.data_1i.copy())
+        self.cube_2i = Cube(self.data_2i.copy())
+        self.cube_1u = Cube(self.data_1u.copy())
+        self.cube_2u = Cube(self.data_2u.copy())
 
         self.ops = (operator.add, operator.sub, operator.mul, operator.truediv)
         self.iops = (
diff --git a/lib/iris/tests/unit/analysis/maths/__init__.py b/lib/iris/tests/unit/analysis/maths/__init__.py
index 521c65a7eb..311da8a0e6 100644
--- a/lib/iris/tests/unit/analysis/maths/__init__.py
+++ b/lib/iris/tests/unit/analysis/maths/__init__.py
@@ -12,6 +12,7 @@
 from abc import ABCMeta, abstractmethod
 import operator
 
+import dask.array as da
 import numpy as np
 from numpy import ma
 
@@ -201,18 +202,22 @@ def cube_func(self):
         # I.E. 'iris.analysis.maths.xx'.
         pass
 
-    def _test_partial_mask(self, in_place):
+    def _test_partial_mask(self, in_place, second_lazy=False):
         # Helper method for masked data tests.
         dat_a = ma.array([2.0, 2.0, 2.0, 2.0], mask=[1, 0, 1, 0])
         dat_b = ma.array([2.0, 2.0, 2.0, 2.0], mask=[1, 1, 0, 0])
 
+        if second_lazy:
+            cube_b = Cube(da.from_array(dat_b))
+        else:
+            cube_b = Cube(dat_b)
+
         cube_a = Cube(dat_a)
-        cube_b = Cube(dat_b)
 
-        com = self.data_op(dat_b, dat_a)
-        res = self.cube_func(cube_b, cube_a, in_place=in_place)
+        com = self.data_op(dat_a, dat_b)
+        res = self.cube_func(cube_a, cube_b, in_place=in_place)
 
-        return com, res, cube_b
+        return com, res, cube_a
 
     def test_partial_mask_in_place(self):
         # Cube in_place arithmetic operation.
@@ -221,13 +226,38 @@ def test_partial_mask_in_place(self):
         self.assertMaskedArrayEqual(com, res.data, strict=True)
         self.assertIs(res, orig_cube)
 
+    def test_partial_mask_second_lazy_in_place(self):
+        # Only second cube has lazy data.
+        com, res, orig_cube = self._test_partial_mask(True, second_lazy=True)
+        self.assertMaskedArrayEqual(com, res.data, strict=True)
+        self.assertIs(res, orig_cube)
+
     def test_partial_mask_not_in_place(self):
         # Cube arithmetic not an in_place operation.
         com, res, orig_cube = self._test_partial_mask(False)
 
-        self.assertMaskedArrayEqual(com, res.data)
+        self.assertMaskedArrayEqual(com, res.data, strict=True)
         self.assertIsNot(res, orig_cube)
 
+    def test_partial_mask_second_lazy_not_in_place(self):
+        # Only second cube has lazy data.
+        com, res, orig_cube = self._test_partial_mask(False, second_lazy=True)
+        self.assertMaskedArrayEqual(com, res.data, strict=True)
+        self.assertIsNot(res, orig_cube)
+
+    def test_in_place_introduces_mask(self):
+        # If second cube is masked, result should also be masked.
+        data1 = np.arange(4, dtype=np.float)
+        data2 = ma.array([2.0, 2.0, 2.0, 2.0], mask=[1, 1, 0, 0])
+        cube1 = Cube(data1)
+        cube2 = Cube(data2)
+
+        com = self.data_op(data1, data2)
+        res = self.cube_func(cube1, cube2, in_place=True)
+
+        self.assertMaskedArrayEqual(com, res.data, strict=True)
+        self.assertIs(res, cube1)
+
 
 class CubeArithmeticCoordsTest(tests.IrisTest):
     # This class sets up pairs of cubes to test iris' ability to reject