diff --git a/CHANGELOG.md b/CHANGELOG.md
index e8f38761..44b797bb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 ### Changed
+- Remove mutable argument defaults to avoid gotchas. ([PR #44](https://github.com/ClimateImpactLab/dscim/pull/44), [@brews](https://github.com/brews))
 - Quiet unused(?), common, logging messages to terminal. ([PR #14](https://github.com/ClimateImpactLab/dscim/pull/14), [@brews](https://github.com/brews))
 
 ## [0.1.0] - 2022-08-30
diff --git a/src/dscim/diagnostics/batch_maps.py b/src/dscim/diagnostics/batch_maps.py
index 6ce065c2..54764b5c 100644
--- a/src/dscim/diagnostics/batch_maps.py
+++ b/src/dscim/diagnostics/batch_maps.py
@@ -74,24 +74,26 @@ def batch_maps(
     ECON_ZARR,
     eta,
     maxmin=True,
-    selection=dict(year=2097, ssp=["SSP2", "SSP3", "SSP4"], batch="batch0"),
+    selection=None,
     location=(0.35, -0.2),
     gcm="mean",
     prefix=None,
-    variables=[
+    variables=(
         "rp_inc_share",
         "rp_damage_share",
         "rp_cons_share",
         "rp",
         "cons",
         "cons_ce",
-    ],
-    maxes=[None, None, None, None, None, None],
-    mins=[None, None, None, None, None, None],
+    ),
+    maxes=(None, None, None, None, None, None),
+    mins=(None, None, None, None, None, None),
     year=2097,
     plot=True,
     save_path=None,
 ):
+    if selection is None:
+        selection = dict(year=2097, ssp=["SSP2", "SSP3", "SSP4"], batch="batch0")
 
     assert len(maxes) == len(
         variables
diff --git a/src/dscim/diagnostics/compare_sccs.py b/src/dscim/diagnostics/compare_sccs.py
index dfb52328..b81f7018 100644
--- a/src/dscim/diagnostics/compare_sccs.py
+++ b/src/dscim/diagnostics/compare_sccs.py
@@ -4,7 +4,7 @@
 
 
 def compare_sccs(
-    rootdict, recipe, disc, eta, rho, quantiles=[0, 0.05, 0.5, 0.95, 1], wp="0.5"
+    rootdict, recipe, disc, eta, rho, quantiles=(0, 0.05, 0.5, 0.95, 1), wp="0.5"
 ):
 
     this_list = []
diff --git a/src/dscim/diagnostics/damage_function.py b/src/dscim/diagnostics/damage_function.py
index c6fa76a3..5c4ccac5 100644
--- a/src/dscim/diagnostics/damage_function.py
+++ b/src/dscim/diagnostics/damage_function.py
@@ -70,7 +70,7 @@ def damage_function(
     rho,
     year=2097,
     hue_vars="ssp",
-    recipes=["adding_up", "risk_aversion", "equity"],
+    recipes=("adding_up", "risk_aversion", "equity"),
     scale=10**12,
     x_lim=(-np.inf, np.inf),
     y_lim=(-np.inf, np.inf),
@@ -97,7 +97,7 @@ def damage_function(
         Type of discounting. Can be 'constant', 'ramsey', or 'wr'.
     year: int
         Year of damage function to be plotted.
-    recipes : list of str
+    recipes : sequence of str
         Recipe types to be plotted. Can be 'adding_up', 'risk_aversion', 'equity'
     scale : int
         Units of dollars for axis. ie., if 10**12 is passed, y axis will be in trillions.
@@ -269,8 +269,8 @@ def damage_function_w_consumption(
     output,
     eta,
     rho,
-    years=[2050, 2099],
-    recipes=["equity"],
+    years=(2050, 2099),
+    recipes=("equity",),
     disc="constant",
     model="IIASA GDP",
     ssp="SSP3",
diff --git a/src/dscim/diagnostics/discount_rates.py b/src/dscim/diagnostics/discount_rates.py
index 9204e849..a58d33af 100644
--- a/src/dscim/diagnostics/discount_rates.py
+++ b/src/dscim/diagnostics/discount_rates.py
@@ -12,13 +12,13 @@ def plot_implicit_rates(
     eta,
     rho,
     fair_aggregation="ce",
-    ssp=["SSP3"],
-    model=["IIASA GDP"],
-    rcp=["ssp585", "ssp245", "ssp460", "ssp370", "Naive Ramsey"],
-    weitzman_parameter=["0.001", "0.1", "1.0"],
+    ssp=("SSP3",),
+    model=("IIASA GDP",),
+    rcp=("ssp585", "ssp245", "ssp460", "ssp370", "Naive Ramsey"),
+    weitzman_parameter=("0.001", "0.1", "1.0"),
     pulse_year=2020,
-    recipes=["risk_aversion", "equity"],
-    discounting=["ramsey", "gwr"],
+    recipes=("risk_aversion", "equity"),
+    discounting=("ramsey", "gwr"),
     aspect=0.4,
     save_path=None,
     csv=True,
diff --git a/src/dscim/diagnostics/equity_risk_premiums.py b/src/dscim/diagnostics/equity_risk_premiums.py
index 178a8b5c..bb3f82d4 100644
--- a/src/dscim/diagnostics/equity_risk_premiums.py
+++ b/src/dscim/diagnostics/equity_risk_premiums.py
@@ -22,18 +22,23 @@ def mumbai_plots(
     eta,
     title=True,
     aspect=1,
-    selection=dict(year=2097),
-    region_dict={
-        "IND.21.317.1249": "Mumbai, IND",
-        "CAN.2.33.913": "Vancouver, CAN",
-        "USA.14.608": "Chicago, USA",
-        "EGY.11": "Cairo, EGY",
-        "SDN.4.11.50.164": "Khartoum, SDN",
-        "NGA.25.510": "Lagos, NGA",
-        "SAU.7": "Riyadh, SAU",
-        "RUS.16.430.430": "St Petersburg, RUS",
-    },
+    selection=None,
+    region_dict=None,
 ):
+    if selection is None:
+        selection = dict(year=2097)
+
+    if region_dict is None:
+        region_dict = {
+            "IND.21.317.1249": "Mumbai, IND",
+            "CAN.2.33.913": "Vancouver, CAN",
+            "USA.14.608": "Chicago, USA",
+            "EGY.11": "Cairo, EGY",
+            "SDN.4.11.50.164": "Khartoum, SDN",
+            "NGA.25.510": "Lagos, NGA",
+            "SAU.7": "Riyadh, SAU",
+            "RUS.16.430.430": "St Petersburg, RUS",
+        }
 
     assert premium in [
         "risk_aversion",
diff --git a/src/dscim/diagnostics/fair_step.py b/src/dscim/diagnostics/fair_step.py
index 8fcfa2ed..a94dde23 100644
--- a/src/dscim/diagnostics/fair_step.py
+++ b/src/dscim/diagnostics/fair_step.py
@@ -21,7 +21,7 @@ def marginal_damages(
     sector,
     sector_path,
     discounting,
-    recipes=["risk_aversion", "equity"],
+    recipes=("risk_aversion", "equity"),
     save_path=None,
 ):
 
@@ -98,7 +98,7 @@ def global_consumption(
     sector,
     sector_path,
     discounting,
-    recipes=["risk_aversion", "equity"],
+    recipes=("risk_aversion", "equity"),
     save_path=None,
     scale=10**12,
 ):
@@ -153,12 +153,12 @@ def output_scc(
     sector_path,
     eta,
     rho,
-    recipes=["adding_up", "risk_aversion", "equity", "local"],
-    discounting=["constant_model_collapsed", "constant", "euler_ramsey", "euler_gwr"],
+    recipes=("adding_up", "risk_aversion", "equity", "local"),
+    discounting=("constant_model_collapsed", "constant", "euler_ramsey", "euler_gwr"),
     save_path=None,
     file=None,
     subset_dict=None,
-    index=[
+    index=(
         "discount_type",
         "discrate",
         "weitzman_parameter",
@@ -166,7 +166,7 @@ def output_scc(
         "ssp",
         "rcp",
         "gas",
-    ],
+    ),
 ):
 
     final_dfs = []
@@ -230,14 +230,19 @@ def plot_implicit_rates(
     sector,
     path,
     fair_aggregation="ce",
-    ssp=["SSP3"],
-    model=["IIASA GDP"],
+    ssp=None,
+    model=None,
     pulse_year=2020,
-    recipes=["risk_aversion", "equity"],
-    discounting=["ramsey", "gwr"],
+    recipes=("risk_aversion", "equity"),
+    discounting=("ramsey", "gwr"),
     save_path=None,
     csv=True,
 ):
+    if ssp is None:
+        ssp = ["SSP3"]
+
+    if model is None:
+        model = ["IIASA GDP"]
 
     for recipe in recipes:
 
diff --git a/src/dscim/diagnostics/maps.py b/src/dscim/diagnostics/maps.py
index 53e2328e..251b031b 100644
--- a/src/dscim/diagnostics/maps.py
+++ b/src/dscim/diagnostics/maps.py
@@ -74,24 +74,26 @@ def maps(
     ECON_ZARR,
     eta,
     maxmin=True,
-    selection=dict(year=2097, ssp=["SSP2", "SSP3", "SSP4"]),
+    selection=None,
     location=(0.35, -0.2),
     gcm="mean",
     prefix=None,
-    variables=[
+    variables=(
         "rp_inc_share",
         "rp_damage_share",
         "rp_cons_share",
         "rp",
         "cons",
         "cons_ce",
-    ],
-    maxes=[None, None, None, None, None, None],
-    mins=[None, None, None, None, None, None],
+    ),
+    maxes=(None, None, None, None, None, None),
+    mins=(None, None, None, None, None, None),
     year=2097,
     plot=True,
     save_path=None,
 ):
+    if selection is None:
+        selection = dict(year=2097, ssp=["SSP2", "SSP3", "SSP4"])
 
     assert len(maxes) == len(
         variables
diff --git a/src/dscim/diagnostics/stacked_damage_function.py b/src/dscim/diagnostics/stacked_damage_function.py
index 545a3374..ab2b176d 100644
--- a/src/dscim/diagnostics/stacked_damage_function.py
+++ b/src/dscim/diagnostics/stacked_damage_function.py
@@ -15,10 +15,12 @@ def plot_stacked(
     eta=2.0,
     rho=0.0,
     xlim=(-1, 8),
-    years=[2020, 2050, 2090, 2100, 2200, 2300],
+    years=None,
     sharey=False,
     rff=True,
 ):
+    if years is None:
+        years = [2020, 2050, 2090, 2100, 2200, 2300]
 
     root_rff = f"/shares/gcp/integration_replication/results/rff/{sector}/2020/"
     root_ssp = f"/shares/gcp/integration_replication/results/AR6_ssp/{sector}/2020/"
diff --git a/src/dscim/menu/main_recipe.py b/src/dscim/menu/main_recipe.py
index f6458a6c..f02f5c9f 100644
--- a/src/dscim/menu/main_recipe.py
+++ b/src/dscim/menu/main_recipe.py
@@ -33,15 +33,18 @@ class MainRecipe(StackedDamages, ABC):
         Discounting is discrete if ``True``, else continuous (default is ``False``).
     fit_type : str
         Type of damage function estimation: ``'ols'``, ``'quantreg'``
-    weitzman_parameter: list of float
+    weitzman_parameter: list of float or None, optional
         If <= 1: The share of global consumption below which bottom coding is implemented.
-        If > 1: Absolute dollar value of global consumption below which bottom
+        If > 1: Absolute dollar value of global consumption below which bottom.
+        Default is [0.1, 0.5].
         coding is implemented.
-    fair_aggregation : list of str
-        How to value climate uncertainty from FAIR: ``median``, ``mean``, ``ce``, ``median_params``
+    fair_aggregation : list of str or None, optional
+        How to value climate uncertainty from FAIR: ``median``, ``mean``,
+        ``ce``, ``median_params``. Default is ["ce", "mean", "gwr_mean",
+        "median", "median_params"].
     rho : float
         Pure rate of time preference parameter
-    fair_dims : list of str
+    fair_dims : list of str or None, optional
         List of dimensions over which the FAIR CE/mean/median options should be collapsed. Default value is ["simulation"], but lists such as ["simulation", "rcp", "ssp"] can be passed. Note: If dimensions other than 'simulation' are passed, 'median_params' fair aggregation cannot be passed.
     """
 
@@ -92,82 +95,105 @@ def __init__(
         damage_function_path=None,
         clip_gmsl=False,
         gdppc_bottom_code=39.39265060424805,
-        scc_quantiles=[0.05, 0.17, 0.25, 0.5, 0.75, 0.83, 0.95],
+        scc_quantiles=None,
         scenario_dimensions=None,
-        weitzman_parameter=[0.1, 0.5],
-        fair_aggregation=["ce", "mean", "gwr_mean", "median", "median_params"],
+        weitzman_parameter=None,
+        fair_aggregation=None,
         filename_suffix="",
         discrete_discounting=False,
-        quantreg_quantiles=[
-            0.05,
-            0.1,
-            0.15,
-            0.2,
-            0.25,
-            0.3,
-            0.35,
-            0.4,
-            0.45,
-            0.5,
-            0.55,
-            0.6,
-            0.65,
-            0.7,
-            0.75,
-            0.8,
-            0.85,
-            0.9,
-            0.95,
-        ],
-        quantreg_weights=[
-            0.075,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.05,
-            0.075,
-        ],
-        full_uncertainty_quantiles=[
-            0.01,
-            0.05,
-            0.17,
-            0.25,
-            0.5,
-            0.75,
-            0.83,
-            0.95,
-            0.99,
-        ],
+        quantreg_quantiles=None,
+        quantreg_weights=None,
+        full_uncertainty_quantiles=None,
         extrap_formula=None,
-        fair_dims=["simulation"],
-        save_files=[
-            "damage_function_points",
-            "damage_function_coefficients",
-            "damage_function_fit",
-            "marginal_damages",
-            "discount_factors",
-            "uncollapsed_sccs",
-            "scc",
-            "uncollapsed_discount_factors",
-            "uncollapsed_marginal_damages",
-            "global_consumption",
-            "global_consumption_no_pulse",
-        ],
+        fair_dims=None,
+        save_files=None,
         **kwargs,
     ):
+        if scc_quantiles is None:
+            scc_quantiles = [0.05, 0.17, 0.25, 0.5, 0.75, 0.83, 0.95]
+
+        if weitzman_parameter is None:
+            weitzman_parameter = [0.1, 0.5]
+
+        if fair_aggregation is None:
+            fair_aggregation = ["ce", "mean", "gwr_mean", "median", "median_params"]
+
+        if quantreg_quantiles is None:
+            quantreg_quantiles = [
+                0.05,
+                0.1,
+                0.15,
+                0.2,
+                0.25,
+                0.3,
+                0.35,
+                0.4,
+                0.45,
+                0.5,
+                0.55,
+                0.6,
+                0.65,
+                0.7,
+                0.75,
+                0.8,
+                0.85,
+                0.9,
+                0.95,
+            ]
+
+        if quantreg_weights is None:
+            quantreg_weights = [
+                0.075,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.05,
+                0.075,
+            ]
+
+        if full_uncertainty_quantiles is None:
+            full_uncertainty_quantiles = [
+                0.01,
+                0.05,
+                0.17,
+                0.25,
+                0.5,
+                0.75,
+                0.83,
+                0.95,
+                0.99,
+            ]
+
+        if fair_dims is None:
+            fair_dims = ["simulation"]
+
+        if save_files is None:
+            save_files = [
+                "damage_function_points",
+                "damage_function_coefficients",
+                "damage_function_fit",
+                "marginal_damages",
+                "discount_factors",
+                "uncollapsed_sccs",
+                "scc",
+                "uncollapsed_discount_factors",
+                "uncollapsed_marginal_damages",
+                "global_consumption",
+                "global_consumption_no_pulse",
+            ]
 
         super().__init__(
             sector_path=sector_path,
diff --git a/src/dscim/menu/simple_storage.py b/src/dscim/menu/simple_storage.py
index a0f03f33..951c0d23 100644
--- a/src/dscim/menu/simple_storage.py
+++ b/src/dscim/menu/simple_storage.py
@@ -33,12 +33,14 @@ class Climate:
         indicating which simulations should be included or excluded.
     ecs_mask_name : str
         Name of mask to be called from within ``ecs_mask_path`` NetCDF file.
-    base_period: tuple
+    base_period: tuple, optional
         Period for rebasing FAIR temperature anomalies. This should match the CIL projection system's base period.
-    emission_scenarios: list
-        List of emission scenarios for which SCC will be calculated.
-    gases: list
-        List of greenhouse gases for which SCC will be calculated.
+    emission_scenarios: list or None, optional
+        List of emission scenarios for which SCC will be calculated. Default
+        is ["ssp119", "ssp126", "ssp245", "ssp460", "ssp370", "ssp585"].
+    gases: list or None, optional
+        List of greenhouse gases for which SCC will be calculated. Default is
+        ["CO2_Fossil", "CH4", "N2O"].
     """
 
     def __init__(
@@ -51,10 +53,22 @@ def __init__(
         gmsl_fair_path=None,
         ecs_mask_path=None,
         ecs_mask_name=None,
-        base_period=[2001, 2010],
-        emission_scenarios=["ssp119", "ssp126", "ssp245", "ssp460", "ssp370", "ssp585"],
-        gases=["CO2_Fossil", "CH4", "N2O"],
+        base_period=(2001, 2010),
+        emission_scenarios=None,
+        gases=None,
     ):
+        if emission_scenarios is None:
+            emission_scenarios = [
+                "ssp119",
+                "ssp126",
+                "ssp245",
+                "ssp460",
+                "ssp370",
+                "ssp585",
+            ]
+        if gases is None:
+            gases = ["CO2_Fossil", "CH4", "N2O"]
+
         self.gmst_path = gmst_path
         self.gmsl_path = gmsl_path
         self.gmst_fair_path = gmst_fair_path
diff --git a/src/dscim/preprocessing/input_damages.py b/src/dscim/preprocessing/input_damages.py
index cae1fa8e..1b5cbcbc 100644
--- a/src/dscim/preprocessing/input_damages.py
+++ b/src/dscim/preprocessing/input_damages.py
@@ -276,7 +276,7 @@ def compute_ag_damages(
     batches=range(0, 15),
     num_cpus=15,
     file="/disaggregated_damages.nc4",
-    vars=["wc_no_reallocation"],
+    vars=None,
     min_year=2010,
     max_year=2099,
 ):
@@ -297,6 +297,8 @@ def compute_ag_damages(
     save_path str
         Path where files should be saved
     """
+    if vars is None:
+        vars = ["wc_no_reallocation"]
 
     if integration:
         assert (
diff --git a/src/dscim/utils/functions.py b/src/dscim/utils/functions.py
index 603920a1..04cea4e7 100644
--- a/src/dscim/utils/functions.py
+++ b/src/dscim/utils/functions.py
@@ -68,7 +68,7 @@ def calculate_constant_equivalent_discount_rate(
     sel_dict,
     eta,
     rho,
-    mean_dims=["simulation"],
+    mean_dims=("simulation",),
     uncollapsed=False,
 ):
     """
diff --git a/src/dscim/utils/menu_runs.py b/src/dscim/utils/menu_runs.py
index 4ae5630c..1a37812f 100644
--- a/src/dscim/utils/menu_runs.py
+++ b/src/dscim/utils/menu_runs.py
@@ -33,10 +33,14 @@ def run_ssps(
     config,
     USA,
     AR,
-    masks=[None],
-    fair_dims_list=[["simulation"]],
+    masks=None,
+    fair_dims_list=None,
     order="damage_function",
 ):
+    if masks is None:
+        masks = [None]
+    if fair_dims_list is None:
+        fair_dims_list = [["simulation"]]
 
     with open(config, "r") as stream:
         conf = yaml.safe_load(stream)
diff --git a/src/dscim/utils/utils.py b/src/dscim/utils/utils.py
index 7408f84d..50cbbc1c 100644
--- a/src/dscim/utils/utils.py
+++ b/src/dscim/utils/utils.py
@@ -103,12 +103,18 @@ def get_weights(quantiles):
     return weights
 
 
-def quantile_weight_quantilereg(
-    array, quantiles=[0.01, 0.05, 0.167, 0.25, 0.5, 0.75, 0.833, 0.95, 0.99]
-):
+def quantile_weight_quantilereg(array, quantiles=None):
     """Produce quantile weights of the quantile regression damages.
+
+    Parameters
+    ----------
     qr_quantiles: the quantile regression quantiles for damages. (Must be 0-1!)
+    quantiles : sequence or None, optional
+        Quantiles to compute. Default is (0.01, 0.05, 0.167, 0.25, 0.5, 0.75,
+        0.833, 0.95, 0.99).
     """
+    if quantiles is None:
+        quantiles = [0.01, 0.05, 0.167, 0.25, 0.5, 0.75, 0.833, 0.95, 0.99]
 
     qr_quantiles = array.q.values
     weights = xr.DataArray(get_weights(qr_quantiles), dims=["q"], coords=[qr_quantiles])
