diff --git a/docs/manual/simulation_configuration.md b/docs/manual/simulation_configuration.md
index f3cc2ccc9ae3bd8632599a00b7d64cf6e8509763..ff4d0a6c581165af51e9d7050cd6bcf0c6877e30 100644
--- a/docs/manual/simulation_configuration.md
+++ b/docs/manual/simulation_configuration.md
@@ -229,6 +229,7 @@ Supported parameters include:
 - Tgrad  
 - Q  
 - fluid_percent  
+- r_out  
 
 Supported sampling distributions and their corresponding arguments for the configuration JSON are:
 
@@ -256,6 +257,10 @@ optimization simulation, see the example for [Optimization of a coaxial BHE](../
 > and upwards in the inner pipe. Optimization of the pipe radius for
 > this design is only supported for the inner pipe.
 
+> When `r_out` is optimized, its resulting value can differ between samples.
+> In that case, `r_out` is stored per-sample (dimensions `samples`; `nPipes`)
+> in the results file instead of the usual single locked-parameter value.
+
 ---
 
 ## Simulation results 
diff --git a/pyproject.toml b/pyproject.toml
index b586f5378bd377cd3773f662f80b133afb4ba369..c653d0b2dda89cb99994bae1da462af0b27aa215 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "geoloop"
-version = "1.0.7"
+version = "1.1.0"
 description = "This is a Python package for simulation of (deep) vertical Borehole Heat Exchanger (BHE) systems"
 authors = [
     { name = "Zanne Korevaar", email = "zanne.korevaar@tno.nl" },
diff --git a/src/geoloop/bin/SingleRunSim.py b/src/geoloop/bin/SingleRunSim.py
index 0be7569db8b72f05084dc39323fef42d0d0f3dab..6f9aa3dfc60948f435cdad3b6f8d8b908da74d53 100644
--- a/src/geoloop/bin/SingleRunSim.py
+++ b/src/geoloop/bin/SingleRunSim.py
@@ -61,7 +61,7 @@ def optimize_forkeys(
         x0[i] = 0.5 * sum(boundsval[i])
 
     method = "Nelder-Mead"
-    method = "Powell"
+    # method = "Powell"
 
     # ftol is tolerance for output power error of optimized point
     result = minimize(
@@ -70,7 +70,7 @@ def optimize_forkeys(
         args=(optimize_keys, config, copcrit, isample),
         bounds=boundsval,
         method=method,
-        options={"maxiter": 50, "ftol": 100},
+        options={"maxiter": 50, "fatol": 5},
     )
 
     return result, config
diff --git a/src/geoloop/constants.py b/src/geoloop/constants.py
index c33539f759bc70ae922204211db7e24a31820c45..4cf84c2a3f4b1692354bcd44e6540d09b737c8ac 100644
--- a/src/geoloop/constants.py
+++ b/src/geoloop/constants.py
@@ -35,6 +35,7 @@ units_dict = {
     "insu_k": "W/mK",
     "insu_dr": "%",
     "k_g_res": "W/mK",
+    "r_out": "m",
 }
 
 # Dictionary of formatted names for various parameters in the model.
@@ -74,6 +75,7 @@ format_dict = {
     "insu_k": "$insu_{k}$",
     "insu_dr": "$insu_{dr}$",
     "k_g_res": "$k_{g-res}$",
+    "r_out": "$r_{out}$",
 }
 
 # Dictionary with colors for lithology plotting
diff --git a/src/geoloop/geoloopcore/CoaxialPipe.py b/src/geoloop/geoloopcore/CoaxialPipe.py
index efec838b6d54e39d8e59fe9629beaed564391bd2..044c35b5b967b9c872b3ba67e4d9afb512bfc33e 100644
--- a/src/geoloop/geoloopcore/CoaxialPipe.py
+++ b/src/geoloop/geoloopcore/CoaxialPipe.py
@@ -198,8 +198,11 @@ class CoaxialPipe(CustomPipe):
         R_f_out_out = 1.0 / (h_f_a_out * 2 * np.pi * r_out_in)
         R_fp = R_p_out + R_f_out_out
 
-        r_inner = np.roll(self.r_in, 1)
-        r_outer = np.roll(self.r_out, 1)
+        # r_inner = np.roll(self.r_in, 1)
+        # r_outer = np.roll(self.r_out, 1)
+
+        r_inner = self.r_in
+        r_outer = self.r_out
 
         self.coaxial = gt.pipes.Coaxial(
             (0, 0), r_inner, r_outer, self.b, k_s, k_g[0], R_ff, R_fp, J=self.J
@@ -498,6 +501,11 @@ class CoaxialPipe(CustomPipe):
         if abs(np.sum(qbz) - power) > 1:
             print("power is not sumq (sumq, power):", np.sum(qbz), ",", power)
 
-        Reff = -np.average(Tf - Tb) * bh.H / power
+        # Reff = -np.average(Tf - Tb) * bh.H / (power + 1e-6)
+
+        if abs(power) > 1e-3:
+            Reff = -np.average(Tf - Tb) * bh.H / power
+        else:
+            Reff = -1
 
         return T_f_out, power, Reff, ptemp, Tb, qbz
diff --git a/src/geoloop/plotting/create_plots.py b/src/geoloop/plotting/create_plots.py
index 937132cafd8dbaa28032304e89545638868c46de..cb966d5346ffc828ff643643ce88b50fcf03af40 100644
--- a/src/geoloop/plotting/create_plots.py
+++ b/src/geoloop/plotting/create_plots.py
@@ -84,16 +84,23 @@ class PlotResults:
         ]
         variable_param_names, locked_param_names = get_param_names()
 
-        x_variables = result_variables + variable_param_names
+        x_variables = result_variables + variable_param_names + ["r_out"]
 
         # Loop through variables and create scatterplots
         for x_variable in x_variables:
             if x_variable != y_variable and x_variable in results_df.columns:
                 # only plot the variables that are truly variable (stdev != 0)
                 if results_df[x_variable].std().round(7) != 0:
-                    title = f"{format_dict[x_variable]} ({units_dict[x_variable]}) after {results_df.loc[0, 'hours']} hours production vs. {format_dict[y_variable]} ({units_dict[y_variable]})"
+                    # Handle r_out unit conversion
+                    plot_df = results_df.copy()
+                    x_unit = units_dict[x_variable]
+                    if x_variable == "r_out":
+                        plot_df[x_variable] = plot_df[x_variable] * 1000
+                        x_unit = "mm"
+
+                    title = f"{format_dict[x_variable]} ({x_unit}) after {plot_df.loc[0, 'hours']} hours production vs. {format_dict[y_variable]} ({units_dict[y_variable]})"
                     g = sns.JointGrid(
-                        results_df, x=x_variable, y=y_data[y_variable], hue="file_name"
+                        plot_df, x=x_variable, y=y_data[y_variable], hue="file_name"
                     )
                     g.plot_joint(
                         sns.scatterplot, alpha=0.7, edgecolor=".2", linewidth=0.5
@@ -118,7 +125,7 @@ class PlotResults:
                         )
                         g.ax_joint.legend()  # Ensure the legend is updated
                     g.set_axis_labels(
-                        xlabel=f"{format_dict[x_variable]} ({units_dict[x_variable]})",
+                        xlabel=f"{format_dict[x_variable]} ({x_unit})",
                         ylabel=f"{format_dict[y_variable]} ({units_dict[y_variable]})",
                     )
                     g.fig.suptitle(title)
@@ -135,9 +142,16 @@ class PlotResults:
                 in params_df.select_dtypes(include=["float64", "int64"]).columns
             ):
                 if params_df[x_variable].std().round(7) != 0:
-                    title = f"{format_dict[x_variable]} ({units_dict[x_variable]}) after {params_df.loc[0, 'nyear']} year production vs. {format_dict[y_variable]} ({units_dict[y_variable]})"
+                    # Handle r_out unit conversion
+                    plot_params_df = params_df.copy()
+                    x_unit = units_dict[x_variable]
+                    if x_variable == "r_out":
+                        plot_params_df[x_variable] = plot_params_df[x_variable] * 1000
+                        x_unit = "mm"
+
+                    title = f"{format_dict[x_variable]} ({x_unit}) after {plot_params_df.loc[0, 'nyear']} year production vs. {format_dict[y_variable]} ({units_dict[y_variable]})"
                     g = sns.JointGrid(
-                        params_df, x=x_variable, y=y_data[y_variable], hue="file_name"
+                        plot_params_df, x=x_variable, y=y_data[y_variable], hue="file_name"
                     )
                     g.plot_joint(
                         sns.scatterplot, alpha=0.7, edgecolor=".2", linewidth=0.5
@@ -156,7 +170,7 @@ class PlotResults:
                         g.ax_joint.invert_yaxis()
 
                     g.set_axis_labels(
-                        xlabel=f"{format_dict[x_variable]} ({units_dict[x_variable]})",
+                        xlabel=f"{format_dict[x_variable]} ({x_unit})",
                         ylabel=f"{format_dict[y_variable]} ({units_dict[y_variable]})",
                     )
                     g.fig.suptitle(title)
diff --git a/src/geoloop/utils/RunManager.py b/src/geoloop/utils/RunManager.py
index d2189918c12ae23433ead269b4f435804b5ccb21..377b228d9056814656d1e85291e2d8e4c2676b21 100644
--- a/src/geoloop/utils/RunManager.py
+++ b/src/geoloop/utils/RunManager.py
@@ -81,10 +81,15 @@ def sample_parameter_space(config: SingleRunConfig) -> tuple[dict, list[str]]:
         sample_space[param_name] = np.full(n_samples, float(fixed_value))
 
     # broadcast locked (non-varying) parameters
+    # Use list() to ensure each sample gets an independent copy of list parameters,
+    # preventing shared-reference issues when optimization modifies values in-place
+    # (e.g. r_out[-1] for COAXIAL pipes)
     for param_name in locked_param_names:
-        sample_space[param_name] = [
-            getattr(config, param_name) for _ in range(n_samples)
-        ]
+        val = getattr(config, param_name)
+        if isinstance(val, list):
+            sample_space[param_name] = [list(val) for _ in range(n_samples)]
+        else:
+            sample_space[param_name] = [val for _ in range(n_samples)]
 
     return sample_space, varying_parameters
 
diff --git a/src/geoloop/utils/helpers.py b/src/geoloop/utils/helpers.py
index 6fc93ec9a73bb3d5cb922c682bae43114b8a6303..3c22f4eaa5604de455b9336931c7fa01cb967d22 100644
--- a/src/geoloop/utils/helpers.py
+++ b/src/geoloop/utils/helpers.py
@@ -625,6 +625,16 @@ def save_MCrun_results(
     variable_param_dict = {key: param_dict[key] for key in variable_param_names}
     locked_param_dict = {key: param_dict[key] for key in locked_param_names}
 
+    # Check if r_out varies across samples (e.g. due to optimization)
+    # by computing the std of the last pipe value across all samples
+    r_out_varies = False
+    if "r_out" in locked_param_dict:
+        r_out_last_values = np.array(
+            [np.asarray(locked_param_dict["r_out"][s])[-1] for s in range(n_samples)]
+        )
+        if round(float(np.std(r_out_last_values)), 7) != 0:
+            r_out_varies = True
+
     # Variable parameters: Create dataset with [param, samples] dimensions
     variable_param_da = xr.DataArray(
         list(variable_param_dict.values()),
@@ -662,16 +672,31 @@ def save_MCrun_results(
     }
 
     pipe_pos = np.asarray(locked_param_x["pos"])
-    r_out = locked_param_x["r_out"]
 
     pipe_pos_ds = xr.Dataset(
         data_vars=dict(pipe_pos=(["nPipes", "xy"], pipe_pos)),
         coords=dict(nPipes=(["nPipes"], pipes_coord)),
     )
-    r_out_ds = xr.Dataset(
-        data_vars=dict(r_out=(["nPipes"], r_out)),
-        coords=dict(nPipes=(["nPipes"], pipes_coord)),
-    )
+
+    # Store r_out: if it varies across samples (optimization), store with [samples, nPipes] dims
+    # otherwise store as locked with [nPipes] dim only
+    if r_out_varies:
+        r_out_per_sample = np.array(
+            [np.asarray(locked_param_dict["r_out"][s]) for s in range(n_samples)]
+        )
+        r_out_ds = xr.Dataset(
+            data_vars=dict(r_out=(["samples", "nPipes"], r_out_per_sample)),
+            coords=dict(samples=(["samples"], list(sample_coord)), nPipes=(["nPipes"], pipes_coord)),
+        )
+    else:
+        r_out = locked_param_x.get("r_out", locked_param_single_dict.get("r_out"))
+        if r_out is not None and isinstance(r_out, list):
+            r_out_ds = xr.Dataset(
+                data_vars=dict(r_out=(["nPipes"], r_out)),
+                coords=dict(nPipes=(["nPipes"], pipes_coord)),
+            )
+        else:
+            r_out_ds = None
 
     # Then, create dataset from locked parameters with dim(len(parnames), nlayers_k_s)
     locked_param_z_k_s = {
diff --git a/uv.lock b/uv.lock
index 4815ea6ac1cb7008a2b16fc3835faa0f2aed8f08..2265c788db49df28ca33d8d26c1a7d6c87fdbbdb 100644
--- a/uv.lock
+++ b/uv.lock
@@ -340,7 +340,7 @@ wheels = [
 
 [[package]]
 name = "geoloop"
-version = "1.0.2"
+version = "1.1.0"
 source = { editable = "." }
 dependencies = [
     { name = "h5netcdf" },