diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index 9c6dff23613549cd7788b05ee5c571af79091179..b8d5c2e3e580820dc210eb6ea0c551082bebbb00 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -6,10 +6,15 @@ from .conversion import bilby_to_lalsimulation_spins
 from .utils import (lalsim_GetApproximantFromString,
                     lalsim_SimInspiralFD,
                     lalsim_SimInspiralChooseFDWaveform,
-                    lalsim_SimInspiralWaveformParamsInsertTidalLambda1,
-                    lalsim_SimInspiralWaveformParamsInsertTidalLambda2,
                     lalsim_SimInspiralChooseFDWaveformSequence)
 
+UNUSED_KWARGS_MESSAGE = """There are unused waveform kwargs. This is deprecated behavior and will
+result in an error in future releases. Make sure all of the waveform kwargs are correctly
+spelled.
+
+Unused waveform_kwargs: {waveform_kwargs}
+"""
+
 
 def gwsignal_binary_black_hole(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
                                phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs):
@@ -480,6 +485,54 @@ def lal_eccentric_binary_black_hole_no_spins(
         eccentricity=eccentricity, **waveform_kwargs)
 
 
+def set_waveform_dictionary(waveform_kwargs, lambda_1=0, lambda_2=0):
+    """
+    Add keyword arguments to the :code:`LALDict` object.
+
+    Parameters
+    ==========
+    waveform_kwargs: dict
+        A dictionary of waveform kwargs. This is modified in place to remove used arguments.
+    lambda_1: float
+        Dimensionless tidal deformability of the primary object.
+    lambda_2: float
+        Dimensionless tidal deformability of the primary object.
+
+    Returns
+    =======
+    waveform_dictionary: lal.LALDict
+        The lal waveform dictionary. This is either taken from the waveform_kwargs or created
+        internally.
+    """
+    import lalsimulation as lalsim
+    from lal import CreateDict
+    waveform_dictionary = waveform_kwargs.pop('lal_waveform_dictionary', CreateDict())
+    waveform_kwargs["TidalLambda1"] = lambda_1
+    waveform_kwargs["TidalLambda2"] = lambda_2
+    waveform_kwargs["NumRelData"] = waveform_kwargs.pop("numerical_relativity_data", None)
+
+    for key in [
+        "pn_spin_order", "pn_tidal_order", "pn_phase_order", "pn_amplitude_order"
+    ]:
+        waveform_kwargs[key[:2].upper() + key[3:].title().replace('_', '')] = waveform_kwargs.pop(key)
+
+    for key in list(waveform_kwargs.keys()).copy():
+        func = getattr(lalsim, f"SimInspiralWaveformParamsInsert{key}", None)
+        if func is None:
+            continue
+        value = waveform_kwargs.pop(key)
+        if func is not None and value is not None:
+            func(waveform_dictionary, value)
+
+    mode_array = waveform_kwargs.pop("mode_array", None)
+    if mode_array is not None:
+        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)
+    return waveform_dictionary
+
+
 def _base_lal_cbc_fd_waveform(
         frequency_array, mass_1, mass_2, luminosity_distance, theta_jn, phase,
         a_1=0.0, a_2=0.0, tilt_1=0.0, tilt_2=0.0, phi_12=0.0, phi_jl=0.0,
@@ -525,22 +578,16 @@ def _base_lal_cbc_fd_waveform(
     =======
     dict: A dictionary with the plus and cross polarisation strain modes
     """
-    import lal
     import lalsimulation as lalsim
 
-    waveform_approximant = waveform_kwargs['waveform_approximant']
-    reference_frequency = waveform_kwargs['reference_frequency']
-    minimum_frequency = waveform_kwargs['minimum_frequency']
-    maximum_frequency = waveform_kwargs['maximum_frequency']
-    catch_waveform_errors = waveform_kwargs['catch_waveform_errors']
-    pn_spin_order = waveform_kwargs['pn_spin_order']
-    pn_tidal_order = waveform_kwargs['pn_tidal_order']
-    pn_phase_order = waveform_kwargs['pn_phase_order']
+    waveform_approximant = waveform_kwargs.pop('waveform_approximant')
+    reference_frequency = waveform_kwargs.pop('reference_frequency')
+    minimum_frequency = waveform_kwargs.pop('minimum_frequency')
+    maximum_frequency = waveform_kwargs.pop('maximum_frequency')
+    catch_waveform_errors = waveform_kwargs.pop('catch_waveform_errors')
     pn_amplitude_order = waveform_kwargs['pn_amplitude_order']
-    waveform_dictionary = waveform_kwargs.get(
-        'lal_waveform_dictionary', lal.CreateDict()
-    )
 
+    waveform_dictionary = set_waveform_dictionary(waveform_kwargs, lambda_1, lambda_2)
     approximant = lalsim_GetApproximantFromString(waveform_approximant)
 
     if pn_amplitude_order != 0:
@@ -567,35 +614,6 @@ def _base_lal_cbc_fd_waveform(
     longitude_ascending_nodes = 0.0
     mean_per_ano = 0.0
 
-    lalsim.SimInspiralWaveformParamsInsertPNSpinOrder(
-        waveform_dictionary, int(pn_spin_order))
-    lalsim.SimInspiralWaveformParamsInsertPNTidalOrder(
-        waveform_dictionary, int(pn_tidal_order))
-    lalsim.SimInspiralWaveformParamsInsertPNPhaseOrder(
-        waveform_dictionary, int(pn_phase_order))
-    lalsim.SimInspiralWaveformParamsInsertPNAmplitudeOrder(
-        waveform_dictionary, int(pn_amplitude_order))
-    lalsim_SimInspiralWaveformParamsInsertTidalLambda1(
-        waveform_dictionary, float(lambda_1))
-    lalsim_SimInspiralWaveformParamsInsertTidalLambda2(
-        waveform_dictionary, float(lambda_2))
-
-    for key, value in waveform_kwargs.items():
-        func = getattr(lalsim, "SimInspiralWaveformParamsInsert" + key, None)
-        if func is not None:
-            func(waveform_dictionary, value)
-
-    if waveform_kwargs.get('numerical_relativity_file', None) is not None:
-        lalsim.SimInspiralWaveformParamsInsertNumRelData(
-            waveform_dictionary, waveform_kwargs['numerical_relativity_file'])
-
-    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
     else:
@@ -650,6 +668,9 @@ def _base_lal_cbc_fd_waveform(
         h_plus[frequency_bounds] *= time_shift
         h_cross[frequency_bounds] *= time_shift
 
+    if len(waveform_kwargs) > 0:
+        logger.warning(UNUSED_KWARGS_MESSAGE.format(waveform_kwargs))
+
     return dict(plus=h_plus, cross=h_cross)
 
 
@@ -705,6 +726,7 @@ def lal_binary_black_hole_relative_binning(
     waveform_kwargs.update(kwargs)
 
     if fiducial == 1:
+        _ = waveform_kwargs.pop("frequency_bin_edges", None)
         return _base_lal_cbc_fd_waveform(
             frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
             luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
@@ -712,6 +734,8 @@ def lal_binary_black_hole_relative_binning(
             phi_12=phi_12, lambda_1=0.0, lambda_2=0.0, **waveform_kwargs)
 
     else:
+        _ = waveform_kwargs.pop("minimum_frequency", None)
+        _ = waveform_kwargs.pop("maximum_frequency", None)
         waveform_kwargs["frequencies"] = waveform_kwargs.pop("frequency_bin_edges")
         return _base_waveform_frequency_sequence(
             frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
@@ -748,6 +772,8 @@ def lal_binary_neutron_star_relative_binning(
             a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_12=phi_12,
             phi_jl=phi_jl, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
     else:
+        _ = waveform_kwargs.pop("minimum_frequency", None)
+        _ = waveform_kwargs.pop("maximum_frequency", None)
         waveform_kwargs["frequencies"] = waveform_kwargs.pop("frequency_bin_edges")
         return _base_waveform_frequency_sequence(
             frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
@@ -822,19 +848,21 @@ def _base_roq_waveform(
     if 'frequency_nodes' not in waveform_arguments:
         size_linear = len(waveform_arguments['frequency_nodes_linear'])
         frequency_nodes_combined = np.hstack(
-            (waveform_arguments['frequency_nodes_linear'],
-             waveform_arguments['frequency_nodes_quadratic'])
+            (waveform_arguments.pop('frequency_nodes_linear'),
+             waveform_arguments.pop('frequency_nodes_quadratic'))
         )
         frequency_nodes_unique, original_indices = np.unique(
             frequency_nodes_combined, return_inverse=True
         )
         linear_indices = original_indices[:size_linear]
         quadratic_indices = original_indices[size_linear:]
-        waveform_arguments['frequency_nodes'] = frequency_nodes_unique
-        waveform_arguments['linear_indices'] = linear_indices
-        waveform_arguments['quadratic_indices'] = quadratic_indices
-
-    waveform_arguments['frequencies'] = waveform_arguments['frequency_nodes']
+        waveform_arguments['frequencies'] = frequency_nodes_unique
+    else:
+        linear_indices = waveform_arguments.pop("linear_indices")
+        quadratic_indices = waveform_arguments.pop("quadratic_indices")
+        for key in ["frequency_nodes_linear", "frequency_nodes_quadratic"]:
+            _ = waveform_arguments.pop(key, None)
+        waveform_arguments['frequencies'] = waveform_arguments.pop('frequency_nodes')
     waveform_polarizations = _base_waveform_frequency_sequence(
         frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
         luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
@@ -843,12 +871,12 @@ def _base_roq_waveform(
 
     return {
         'linear': {
-            'plus': waveform_polarizations['plus'][waveform_arguments['linear_indices']],
-            'cross': waveform_polarizations['cross'][waveform_arguments['linear_indices']]
+            'plus': waveform_polarizations['plus'][linear_indices],
+            'cross': waveform_polarizations['cross'][linear_indices]
         },
         'quadratic': {
-            'plus': waveform_polarizations['plus'][waveform_arguments['quadratic_indices']],
-            'cross': waveform_polarizations['cross'][waveform_arguments['quadratic_indices']]
+            'plus': waveform_polarizations['plus'][quadratic_indices],
+            'cross': waveform_polarizations['cross'][quadratic_indices]
         }
     }
 
@@ -1059,49 +1087,13 @@ def _base_waveform_frequency_sequence(
         Dict containing plus and cross modes evaluated at the linear and
         quadratic frequency nodes.
     """
-    from lal import CreateDict
-    import lalsimulation as lalsim
-
-    frequencies = waveform_kwargs['frequencies']
-    reference_frequency = waveform_kwargs['reference_frequency']
-    approximant = lalsim_GetApproximantFromString(waveform_kwargs['waveform_approximant'])
-    catch_waveform_errors = waveform_kwargs['catch_waveform_errors']
-    pn_spin_order = waveform_kwargs['pn_spin_order']
-    pn_tidal_order = waveform_kwargs['pn_tidal_order']
-    pn_phase_order = waveform_kwargs['pn_phase_order']
-    pn_amplitude_order = waveform_kwargs['pn_amplitude_order']
-    waveform_dictionary = waveform_kwargs.get(
-        'lal_waveform_dictionary', CreateDict()
-    )
-
-    lalsim.SimInspiralWaveformParamsInsertPNSpinOrder(
-        waveform_dictionary, int(pn_spin_order))
-    lalsim.SimInspiralWaveformParamsInsertPNTidalOrder(
-        waveform_dictionary, int(pn_tidal_order))
-    lalsim.SimInspiralWaveformParamsInsertPNPhaseOrder(
-        waveform_dictionary, int(pn_phase_order))
-    lalsim.SimInspiralWaveformParamsInsertPNAmplitudeOrder(
-        waveform_dictionary, int(pn_amplitude_order))
-    lalsim_SimInspiralWaveformParamsInsertTidalLambda1(
-        waveform_dictionary, float(lambda_1))
-    lalsim_SimInspiralWaveformParamsInsertTidalLambda2(
-        waveform_dictionary, float(lambda_2))
-
-    for key, value in waveform_kwargs.items():
-        func = getattr(lalsim, "SimInspiralWaveformParamsInsert" + key, None)
-        if func is not None:
-            func(waveform_dictionary, value)
+    frequencies = waveform_kwargs.pop('frequencies')
+    reference_frequency = waveform_kwargs.pop('reference_frequency')
+    approximant = waveform_kwargs.pop('waveform_approximant')
+    catch_waveform_errors = waveform_kwargs.pop('catch_waveform_errors')
 
-    if waveform_kwargs.get('numerical_relativity_file', None) is not None:
-        lalsim.SimInspiralWaveformParamsInsertNumRelData(
-            waveform_dictionary, waveform_kwargs['numerical_relativity_file'])
-
-    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)
+    waveform_dictionary = set_waveform_dictionary(waveform_kwargs, lambda_1, lambda_2)
+    approximant = lalsim_GetApproximantFromString(approximant)
 
     luminosity_distance = luminosity_distance * 1e6 * utils.parsec
     mass_1 = mass_1 * utils.solar_mass
@@ -1135,6 +1127,9 @@ def _base_waveform_frequency_sequence(
             else:
                 raise
 
+    if len(waveform_kwargs) > 0:
+        logger.warning(UNUSED_KWARGS_MESSAGE.format(waveform_kwargs))
+
     return dict(plus=h_plus.data.data, cross=h_cross.data.data)
 
 
diff --git a/test/gw/likelihood/marginalization_test.py b/test/gw/likelihood/marginalization_test.py
index d7b4f8a59785e695b5b6d58cb5f036c60b4e5101..498ca9fe1fb9f9f6b7fd42956a0a6dc192efd0f6 100644
--- a/test/gw/likelihood/marginalization_test.py
+++ b/test/gw/likelihood/marginalization_test.py
@@ -248,7 +248,6 @@ class TestMarginalizations(unittest.TestCase):
             start_time=1126259640,
             waveform_arguments=dict(
                 reference_frequency=20.0,
-                minimum_frequency=20.0,
                 waveform_approximant="IMRPhenomPv2",
                 frequency_nodes_linear=np.load(f"{roq_dir}/fnodes_linear.npy"),
                 frequency_nodes_quadratic=np.load(f"{roq_dir}/fnodes_quadratic.npy"),
diff --git a/test/gw/likelihood_test.py b/test/gw/likelihood_test.py
index b12ffd59bdc9101bca7170c5764cd73c7a9eecbb..84298e1d9f6a7cf48b7bde1bb42d50ca4b49790b 100644
--- a/test/gw/likelihood_test.py
+++ b/test/gw/likelihood_test.py
@@ -361,7 +361,6 @@ class TestROQLikelihood(unittest.TestCase):
                 frequency_nodes_linear=fnodes_linear,
                 frequency_nodes_quadratic=fnodes_quadratic,
                 reference_frequency=20.0,
-                minimum_frequency=20.0,
                 waveform_approximant="IMRPhenomPv2",
             ),
         )
@@ -599,7 +598,6 @@ class TestRescaledROQLikelihood(unittest.TestCase):
                 frequency_nodes_linear=fnodes_linear,
                 frequency_nodes_quadratic=fnodes_quadratic,
                 reference_frequency=20.0,
-                minimum_frequency=20.0,
                 waveform_approximant="IMRPhenomPv2",
             ),
         )
diff --git a/test/gw/source_test.py b/test/gw/source_test.py
index 5f5199fd10698993036f21fdfe5517cbf7dcbff7..8f55c5b14be1151a456ad4dcf5f751346f5b8ec5 100644
--- a/test/gw/source_test.py
+++ b/test/gw/source_test.py
@@ -305,7 +305,6 @@ class TestROQBBH(unittest.TestCase):
             frequency_nodes_linear=fnodes_linear,
             frequency_nodes_quadratic=fnodes_quadratic,
             reference_frequency=50.0,
-            minimum_frequency=20.0,
             waveform_approximant="IMRPhenomPv2",
         )
         self.frequency_array = bilby.core.utils.create_frequency_series(2048, 4)
@@ -344,15 +343,14 @@ class TestBBHfreqseq(unittest.TestCase):
             theta_jn=0.3,
             phase=0.0
         )
-        minimum_frequency = 20.0
+        self.minimum_frequency = 20.0
         self.frequency_array = bilby.core.utils.create_frequency_series(2048, 8)
-        self.full_frequencies_to_sequence = self.frequency_array >= minimum_frequency
+        self.full_frequencies_to_sequence = self.frequency_array >= self.minimum_frequency
+        self.frequencies = self.frequency_array[self.full_frequencies_to_sequence]
         self.waveform_kwargs = dict(
             waveform_approximant="IMRPhenomHM",
             reference_frequency=50.0,
-            minimum_frequency=minimum_frequency,
             catch_waveform_errors=True,
-            frequencies=self.frequency_array[self.full_frequencies_to_sequence]
         )
         self.bad_parameters = copy(self.parameters)
         self.bad_parameters["mass_1"] = -30.0
@@ -362,12 +360,13 @@ class TestBBHfreqseq(unittest.TestCase):
         del self.waveform_kwargs
         del self.frequency_array
         del self.bad_parameters
+        del self.minimum_frequency
 
     def test_valid_parameters(self):
         self.parameters.update(self.waveform_kwargs)
         self.assertIsInstance(
             bilby.gw.source.binary_black_hole_frequency_sequence(
-                self.frequency_array, **self.parameters
+                self.frequency_array, frequencies=self.frequencies, **self.parameters
             ),
             dict
         )
@@ -376,7 +375,7 @@ class TestBBHfreqseq(unittest.TestCase):
         self.bad_parameters.update(self.waveform_kwargs)
         self.assertIsNone(
             bilby.gw.source.binary_black_hole_frequency_sequence(
-                self.frequency_array, **self.bad_parameters
+                self.frequency_array, frequencies=self.frequencies, **self.bad_parameters
             )
         )
 
@@ -386,16 +385,16 @@ class TestBBHfreqseq(unittest.TestCase):
         raise_error_parameters["catch_waveform_errors"] = False
         with self.assertRaises(Exception):
             bilby.gw.source.binary_black_hole_frequency_sequence(
-                self.frequency_array, **raise_error_parameters
+                self.frequency_array, frequencies=self.frequencies, **raise_error_parameters
             )
 
     def test_match_LalBBH(self):
         self.parameters.update(self.waveform_kwargs)
         freqseq = bilby.gw.source.binary_black_hole_frequency_sequence(
-            self.frequency_array, **self.parameters
+            self.frequency_array, frequencies=self.frequencies, **self.parameters
         )
         lalbbh = bilby.gw.source.lal_binary_black_hole(
-            self.frequency_array, **self.parameters
+            self.frequency_array, minimum_frequency=self.minimum_frequency, **self.parameters
         )
         self.assertEqual(freqseq.keys(), lalbbh.keys())
         for mode in freqseq:
@@ -408,10 +407,10 @@ class TestBBHfreqseq(unittest.TestCase):
         parameters.update(self.waveform_kwargs)
         parameters['mode_array'] = [[2, 2]]
         freqseq = bilby.gw.source.binary_black_hole_frequency_sequence(
-            self.frequency_array, **parameters
+            self.frequency_array, frequencies=self.frequencies, **parameters
         )
         lalbbh = bilby.gw.source.lal_binary_black_hole(
-            self.frequency_array, **parameters
+            self.frequency_array, minimum_frequency=self.minimum_frequency, **parameters
         )
         self.assertEqual(freqseq.keys(), lalbbh.keys())
         for mode in freqseq:
@@ -426,10 +425,10 @@ class TestBBHfreqseq(unittest.TestCase):
         lalsimulation.SimInspiralWaveformParamsInsertNonGRDChi0(wf_dict, 1.)
         parameters['lal_waveform_dictionary'] = wf_dict
         freqseq = bilby.gw.source.binary_black_hole_frequency_sequence(
-            self.frequency_array, **parameters
+            self.frequency_array, frequencies=self.frequencies, **parameters
         )
         lalbbh = bilby.gw.source.lal_binary_black_hole(
-            self.frequency_array, **parameters
+            self.frequency_array, minimum_frequency=self.minimum_frequency, **parameters
         )
         self.assertEqual(freqseq.keys(), lalbbh.keys())
         for mode in freqseq:
@@ -455,26 +454,26 @@ class TestBNSfreqseq(unittest.TestCase):
             lambda_1=1000.0,
             lambda_2=1000.0
         )
-        minimum_frequency = 50.0
+        self.minimum_frequency = 50.0
         self.frequency_array = bilby.core.utils.create_frequency_series(2048, 16)
-        self.full_frequencies_to_sequence = self.frequency_array >= minimum_frequency
+        self.full_frequencies_to_sequence = self.frequency_array >= self.minimum_frequency
+        self.frequencies = self.frequency_array[self.full_frequencies_to_sequence]
         self.waveform_kwargs = dict(
             waveform_approximant="IMRPhenomPv2_NRTidal",
             reference_frequency=50.0,
-            minimum_frequency=minimum_frequency,
-            frequencies=self.frequency_array[self.full_frequencies_to_sequence]
         )
 
     def tearDown(self):
         del self.parameters
         del self.waveform_kwargs
         del self.frequency_array
+        del self.minimum_frequency
 
     def test_with_valid_parameters(self):
         self.parameters.update(self.waveform_kwargs)
         self.assertIsInstance(
             bilby.gw.source.binary_neutron_star_frequency_sequence(
-                self.frequency_array, **self.parameters
+                self.frequency_array, frequencies=self.frequencies, **self.parameters
             ),
             dict
         )
@@ -485,16 +484,16 @@ class TestBNSfreqseq(unittest.TestCase):
         self.parameters.update(self.waveform_kwargs)
         with self.assertRaises(TypeError):
             bilby.gw.source.binary_neutron_star_frequency_sequence(
-                self.frequency_array, **self.parameters
+                self.frequency_array, frequencies=self.frequencies, **self.parameters
             )
 
     def test_match_LalBNS(self):
         self.parameters.update(self.waveform_kwargs)
         freqseq = bilby.gw.source.binary_neutron_star_frequency_sequence(
-            self.frequency_array, **self.parameters
+            self.frequency_array, frequencies=self.frequencies, **self.parameters
         )
         lalbns = bilby.gw.source.lal_binary_neutron_star(
-            self.frequency_array, **self.parameters
+            self.frequency_array, minimum_frequency=self.minimum_frequency, **self.parameters
         )
         self.assertEqual(freqseq.keys(), lalbns.keys())
         for mode in freqseq: