diff --git a/src/bayesvalidrox/post_processing/post_processing.py b/src/bayesvalidrox/post_processing/post_processing.py
index bed9851055ccd41df8d67fe17e5d99d96ab9b3f8..b29d8bf40b06dd9a76c92912ad92092407e7a7c5 100644
--- a/src/bayesvalidrox/post_processing/post_processing.py
+++ b/src/bayesvalidrox/post_processing/post_processing.py
@@ -230,7 +230,7 @@ class PostProcessing:
             Parameter sets to be checked. The default is None.
         outputs : dict, optional
             Output dictionary with model outputs for all given output types in
-            `Model.Output.names`. The default is None.
+            `engine.out_names`. The default is None.
 
         Raises
         ------
@@ -555,7 +555,7 @@ class PostProcessing:
 
         # Extract the necessary variables
         max_order = np.max(metamod._pce_deg)
-        outputs = self.engine.Model.Output.names
+        outputs = self.engine.out_names
         if metamod.sobol is None:
             metamod.calculate_sobol()
         sobol_all, total_sobol_all = metamod.sobol, metamod.total_sobol
@@ -700,6 +700,7 @@ class PostProcessing:
             )
             plt.clf()
 
+        return sobol_all, total_sobol_all
 
     # -------------------------------------------------------------------------
     def check_reg_quality(self, n_samples: int = 1000, samples=None, outputs: dict = None) -> None:
@@ -715,7 +716,7 @@ class PostProcessing:
             Parameter sets to use for the check. The default is None.
         outputs : dict, optional
             Output dictionary with model outputs for all given output types in
-            `Model.Output.names`. The default is None.
+            `engine.out_names`. The default is None.
 
         Return 
         ------
diff --git a/src/bayesvalidrox/surrogate_models/polynomial_chaos.py b/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
index afedbe76440300104ae0f4572d8f47e35ad1071e..c1164b690eb292a536abde52823e9ec90ae56fdc 100644
--- a/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
+++ b/src/bayesvalidrox/surrogate_models/polynomial_chaos.py
@@ -929,7 +929,7 @@ class PCE(MetaModel):
 
         return pce_means, pce_stds
 
-    def calculate_sobol(self, Y_train = None):
+    def calculate_sobol(self, y_train = None):
         """
         Provides Sobol' indices as a sensitivity measure to infer the importance
         of the input parameters. See Eq. 27 in [1] for more details. For the
@@ -947,7 +947,7 @@ class PCE(MetaModel):
 
         Parameters
         ----------
-        Y_train: dict, optional
+        y_train: dict, optional
             Trainings outputs. They are needed when used in combination with PCA.
             The default is None
 
@@ -959,6 +959,9 @@ class PCE(MetaModel):
             Total Sobol' indices
 
         """
+        if self.dim_red_method == 'pca' and y_train is None:
+            raise AttributeError("Calculation of Sobol' indices with PCA expects training outputs, but none are given.")
+
         max_order = np.max(self._pce_deg)
         n_params = len(self.input_obj.Marginals)
         cov_zpq = np.zeros((n_params))
@@ -1061,7 +1064,7 @@ class PCE(MetaModel):
 
                 # Compute the sobol indices according to Ref. 2
                 if self.dim_red_method.lower() == "pca":
-                    n_c_points = Y_train[output].shape[1]
+                    n_c_points = y_train[output].shape[1]
                     pca = self.pca[f"b_{b_i+1}"][output]
                     comp_pca = pca.components_
                     n_comp = comp_pca.shape[0]
@@ -1079,7 +1082,7 @@ class PCE(MetaModel):
 
                             for t_idx in range(n_c_points):
                                 var_yt = np.var(
-                                    Y_train[output][:, t_idx]
+                                    y_train[output][:, t_idx]
                                 )
                                 term1, term2 = 0.0, 0.0
                                 if var_yt != 0.0:
@@ -1106,7 +1109,7 @@ class PCE(MetaModel):
                         s_zi = total_sobol_array[par_idx]
 
                         for t_idx in range(n_c_points):
-                            var_yt = np.var(Y_train[output][:, t_idx])
+                            var_yt = np.var(y_train[output][:, t_idx])
                             term1, term2 = 0.0, 0.0
                             if var_yt != 0.0:
                                 for i in range(n_comp):
diff --git a/tests/test_PolynomialChaosEmulator.py b/tests/test_PolynomialChaosEmulator.py
index 64858588d1382caa8df9002fec2ed4c2d7ef9b9a..59477eda68810a68e6aca8c9bb7ce855ce2f674d 100644
--- a/tests/test_PolynomialChaosEmulator.py
+++ b/tests/test_PolynomialChaosEmulator.py
@@ -873,3 +873,20 @@ def test_calculate_sobol_pca():
     mm.fit([[0.2], [0.8]], {'Z': [[0.4, 0.4], [0.5, 0.6]]})
     sobol, totalsobol = mm.calculate_sobol({'Z': np.array([[0.4, 0.4], [0.5, 0.6]])})
     # TODO are there theory-related checks that could be applied here?
+
+def test_calculate_sobol_pcanoy():
+    """
+    Calculate Sobol' indices of a pce-surrogate with PCA but no outputs
+    """
+    inp = Input()
+    inp.add_marginals()
+    inp.Marginals[0].dist_type = 'normal'
+    inp.Marginals[0].parameters = [0, 1]
+    mm = PCE(inp)
+    mm.dim_red_method = 'pca'
+    mm.fit([[0.2], [0.8]], {'Z': [[0.4, 0.4], [0.5, 0.6]]})
+
+    with pytest.raises(AttributeError) as excinfo:
+        mm.calculate_sobol()
+    assert str(excinfo.value) == ("Calculation of Sobol' indices with PCA expects training outputs, but none are given.")
+
diff --git a/tests/test_PostProcessing.py b/tests/test_PostProcessing.py
index a57e8c085cc3f58b20771768743de87ffbec7a24..eb3d194663a3910eba545ee13051a3227180d65b 100644
--- a/tests/test_PostProcessing.py
+++ b/tests/test_PostProcessing.py
@@ -143,10 +143,17 @@ def test_sobol_indices_pce(pce_engine) -> None:
     """
     engine = pce_engine
     post = PostProcessing(engine)
-    sobol = post.sobol_indices()
-    assert list(sobol.keys()) == ['Z']
-    assert sobol['Z'].shape == (1,1)
-    assert sobol['Z'][0,0] == 1
+    sobol, totalsobol = post.sobol_indices()
+
+    assert list(totalsobol.keys()) == ['Z']
+    assert totalsobol['Z'].shape == (1,1)
+    assert totalsobol['Z'][0,0] == 1
+
+    print(sobol)
+    assert list(sobol.keys()) == [1]
+    assert list(sobol[1].keys()) == ['Z']
+    assert sobol[1]['Z'].shape == (1,1,1)
+    assert sobol[1]['Z'][0,0] == 1
 
 #%% check_reg_quality