diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index 31271f8a810371564eb8c13569e92b96fe431547..69109364a149304926318ebd8b73fc1adf96425e 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -67,18 +67,18 @@ def lal_binary_black_hole(
             mode_array:
                 Activate a specific mode array and evaluate the model using those
                 modes only.  e.g. waveform_arguments =
-                dict(waveform_approximant='IMRPhenomHM', modearray=[[2,2],[2,-2])
+                dict(waveform_approximant='IMRPhenomHM', mode_array=[[2,2],[2,-2])
                 returns the 22 and 2-2 modes only of IMRPhenomHM.  You can only
                 specify modes that are included in that particular model.  e.g.
                 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
-                modearray=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
+                mode_array=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
                 55 modes are not included in this model.  Be aware that some models
                 only take positive modes and return the positive and the negative
                 mode together, while others need to call both.  e.g.
                 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
-                modearray=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
+                mode_array=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
                 However, waveform_arguments =
-                dict(waveform_approximant='IMRPhenomXHM', modearray=[[2,2],[4,-4]])
+                dict(waveform_approximant='IMRPhenomXHM', mode_array=[[2,2],[4,-4]])
                 returns the 22 and 4-4 of IMRPhenomXHM.
 
     Returns
@@ -150,18 +150,18 @@ def lal_binary_neutron_star(
             mode_array:
                 Activate a specific mode array and evaluate the model using those
                 modes only.  e.g. waveform_arguments =
-                dict(waveform_approximant='IMRPhenomHM', modearray=[[2,2],[2,-2])
+                dict(waveform_approximant='IMRPhenomHM', mode_array=[[2,2],[2,-2])
                 returns the 22 and 2-2 modes only of IMRPhenomHM.  You can only
                 specify modes that are included in that particular model.  e.g.
                 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
-                modearray=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
+                mode_array=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
                 55 modes are not included in this model.  Be aware that some models
                 only take positive modes and return the positive and the negative
                 mode together, while others need to call both.  e.g.
                 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
-                modearray=[[2,2],[4,-4]]) returns the 22 a\nd 2-2 of IMRPhenomHM.
+                mode_array=[[2,2],[4,-4]]) returns the 22 a\nd 2-2 of IMRPhenomHM.
                 However, waveform_arguments =
-                dict(waveform_approximant='IMRPhenomXHM', modearray=[[2,2],[4,-4]])
+                dict(waveform_approximant='IMRPhenomXHM', mode_array=[[2,2],[4,-4]])
                 returns the 22 and 4-4 of IMRPhenomXHM.
 
     Returns
@@ -217,18 +217,18 @@ def lal_eccentric_binary_black_hole_no_spins(
             mode_array:
                 Activate a specific mode array and evaluate the model using those
                 modes only.  e.g. waveform_arguments =
-                dict(waveform_approximant='IMRPhenomHM', modearray=[[2,2],[2,-2])
+                dict(waveform_approximant='IMRPhenomHM', mode_array=[[2,2],[2,-2])
                 returns the 22 and 2-2 modes only of IMRPhenomHM.  You can only
                 specify modes that are included in that particular model.  e.g.
                 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
-                modearray=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
+                mode_array=[[2,2],[2,-2],[5,5],[5,-5]]) is not allowed because the
                 55 modes are not included in this model.  Be aware that some models
                 only take positive modes and return the positive and the negative
                 mode together, while others need to call both.  e.g.
                 waveform_arguments = dict(waveform_approximant='IMRPhenomHM',
-                modearray=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
+                mode_array=[[2,2],[4,-4]]) returns the 22 and 2-2 of IMRPhenomHM.
                 However, waveform_arguments =
-                dict(waveform_approximant='IMRPhenomXHM', modearray=[[2,2],[4,-4]])
+                dict(waveform_approximant='IMRPhenomXHM', mode_array=[[2,2],[4,-4]])
                 returns the 22 and 4-4 of IMRPhenomXHM.
 
     Returns
@@ -343,12 +343,12 @@ def _base_lal_cbc_fd_waveform(
     lalsim_SimInspiralWaveformParamsInsertTidalLambda2(
         waveform_dictionary, lambda_2)
 
-    if 'modearray' in waveform_kwargs:
-        modearray = waveform_kwargs['modearray']
-        mode_array = lalsim.SimInspiralCreateModeArray()
-        for mode in modearray:
-            lalsim.SimInspiralModeArrayActivateMode(mode_array, mode[0], mode[1])
-        lalsim.SimInspiralWaveformParamsInsertModeArray(waveform_dictionary, mode_array)
+    if ('mode_array' in waveform_kwargs) and waveform_kwargs['mode_array'] is not None:
+        mode_array = waveform_kwargs['mode_array']
+        mode_array_lal = lalsim.SimInspiralCreateModeArray()
+        for mode in mode_array:
+            lalsim.SimInspiralModeArrayActivateMode(mode_array_lal, mode[0], mode[1])
+        lalsim.SimInspiralWaveformParamsInsertModeArray(waveform_dictionary, mode_array_lal)
 
     if lalsim.SimInspiralImplementedFDApproximants(approximant):
         wf_func = lalsim_SimInspiralChooseFDWaveform