diff --git a/include/amici/defines.h b/include/amici/defines.h
index b49ebc6a83..44fa4cadc5 100644
--- a/include/amici/defines.h
+++ b/include/amici/defines.h
@@ -244,6 +244,7 @@ enum class RDataReporting {
     full,
     residuals,
     likelihood,
+    observables_likelihood,
 };
 
 /** boundary conditions for splines */
diff --git a/include/amici/rdata.h b/include/amici/rdata.h
index 793be9435a..9f7378751a 100644
--- a/include/amici/rdata.h
+++ b/include/amici/rdata.h
@@ -476,6 +476,13 @@ class ReturnData : public ModelDimensions {
      */
     void initializeLikelihoodReporting(bool quadratic_llh);
 
+    /**
+     * @brief initializes storage for observables + likelihood reporting mode
+     * @param quadratic_llh whether model defines a quadratic nllh and computing
+     * res, sres and FIM makes sense.
+     */
+    void initializeObservablesLikelihoodReporting(bool quadratic_llh);
+
     /**
      * @brief initializes storage for residual reporting mode
      * @param enable_res whether residuals are to be computed
diff --git a/python/tests/test_swig_interface.py b/python/tests/test_swig_interface.py
index 9686a25d94..4ead568c9c 100644
--- a/python/tests/test_swig_interface.py
+++ b/python/tests/test_swig_interface.py
@@ -586,3 +586,44 @@ def test_python_exceptions(sbml_example_presimulation_module):
     ):
         # rethrow=True
         runAmiciSimulation(solver, None, model.get(), True)
+
+
+def test_reporting_mode_obs_llh(sbml_example_presimulation_module):
+    model_module = sbml_example_presimulation_module
+    model = model_module.getModel()
+    solver = model.getSolver()
+
+    solver.setReturnDataReportingMode(
+        amici.RDataReporting.observables_likelihood
+    )
+    solver.setSensitivityOrder(amici.SensitivityOrder.first)
+
+    for sens_method in (
+        amici.SensitivityMethod.none,
+        amici.SensitivityMethod.forward,
+        amici.SensitivityMethod.adjoint,
+    ):
+        solver.setSensitivityMethod(sens_method)
+        rdata = amici.runAmiciSimulation(
+            model, solver, amici.ExpData(1, 1, 1, [1])
+        )
+        assert (
+            rdata.rdata_reporting
+            == amici.RDataReporting.observables_likelihood
+        )
+
+        assert rdata.y.size > 0
+        assert rdata.sigmay.size > 0
+        assert rdata.J is None
+
+        match solver.getSensitivityMethod():
+            case amici.SensitivityMethod.none:
+                assert rdata.sllh is None
+            case amici.SensitivityMethod.forward:
+                assert rdata.sy.size > 0
+                assert rdata.ssigmay.size > 0
+                assert rdata.sllh.size > 0
+            case amici.SensitivityMethod.adjoint:
+                assert rdata.sy is None
+                assert rdata.ssigmay is None
+                assert rdata.sllh.size > 0
diff --git a/src/rdata.cpp b/src/rdata.cpp
index 4ec983af2b..c724d29954 100644
--- a/src/rdata.cpp
+++ b/src/rdata.cpp
@@ -60,13 +60,17 @@ ReturnData::ReturnData(
     case RDataReporting::likelihood:
         initializeLikelihoodReporting(quadratic_llh);
         break;
+
+    case RDataReporting::observables_likelihood:
+        initializeObservablesLikelihoodReporting(quadratic_llh);
+        break;
     }
 }
 
 void ReturnData::initializeLikelihoodReporting(bool enable_fim) {
     llh = getNaN();
     chi2 = getNaN();
-    if (sensi >= SensitivityOrder::first) {
+    if (sensi >= SensitivityOrder::first && sensi_meth != SensitivityMethod::none) {
         sllh.resize(nplist, getNaN());
         if (sensi >= SensitivityOrder::second)
             s2llh.resize(nplist * (nJ - 1), getNaN());
@@ -78,6 +82,21 @@ void ReturnData::initializeLikelihoodReporting(bool enable_fim) {
     }
 }
 
+void ReturnData::initializeObservablesLikelihoodReporting(bool enable_fim) {
+    initializeLikelihoodReporting(enable_fim);
+
+    y.resize(nt * ny, 0.0);
+    sigmay.resize(nt * ny, 0.0);
+
+    if ((sensi_meth == SensitivityMethod::forward
+         && sensi >= SensitivityOrder::first)
+        || sensi >= SensitivityOrder::second) {
+
+        sy.resize(nt * ny * nplist, 0.0);
+        ssigmay.resize(nt * ny * nplist, 0.0);
+    }
+}
+
 void ReturnData::initializeResidualReporting(bool enable_res) {
     y.resize(nt * ny, 0.0);
     sigmay.resize(nt * ny, 0.0);