diff --git a/src/storm-cli-utilities/model-handling.h b/src/storm-cli-utilities/model-handling.h
index cc7c875c4..efeb46345 100644
--- a/src/storm-cli-utilities/model-handling.h
+++ b/src/storm-cli-utilities/model-handling.h
@@ -702,7 +702,8 @@ std::shared_ptr<storm::models::sparse::Model<ValueType>> preprocessSparseModelBi
     std::optional<double> tolerance = storm::settings::getModule<storm::settings::modules::GeneralSettings>().getPrecision();
 
     STORM_LOG_INFO("Performing bisimulation minimization...");
-    return storm::api::performBisimulationMinimization<ValueType>(model, createFormulasToRespect(input.properties), bisimType, graphPreserving, tolerance);
+    return storm::api::performBisimulationMinimization<ValueType>(model, createFormulasToRespect(input.properties), bisimType, graphPreserving, tolerance,
+                                                                  !bisimulationSettings.isMeasureDrivenPartitionDisabled());
 }
 
 template<typename ValueType>
diff --git a/src/storm/api/bisimulation.h b/src/storm/api/bisimulation.h
index 074cfb32d..b34780404 100644
--- a/src/storm/api/bisimulation.h
+++ b/src/storm/api/bisimulation.h
@@ -24,7 +24,8 @@ template<typename ModelType>
 std::shared_ptr<ModelType> performDeterministicSparseBisimulationMinimization(std::shared_ptr<ModelType> model,
                                                                               std::vector<std::shared_ptr<storm::logic::Formula const>> const& formulas,
                                                                               storm::storage::BisimulationType type, bool graphPreserving = true,
-                                                                              std::optional<double> const& tolerance = std::nullopt) {
+                                                                              std::optional<double> const& tolerance = std::nullopt,
+                                                                              bool allowMeasureDrivenInitialPartition = true) {
     using OptionsType = typename storm::storage::DeterministicModelBisimulationDecomposition<ModelType>::Options;
     // Falls back to the general precision setting when the caller does not deliberately choose a tolerance;
     // may be reworked to require an explicit choice throughout the API in the future.
@@ -42,6 +43,7 @@ std::shared_ptr<ModelType> performDeterministicSparseBisimulationMinimization(st
         options.setKeepRewards(true);
     }
     options.setType(type);
+    options.setAllowMeasureDrivenInitialPartition(allowMeasureDrivenInitialPartition);
 
     storm::storage::DeterministicModelBisimulationDecomposition<ModelType> bisimulationDecomposition(*model, options);
     bisimulationDecomposition.computeBisimulationDecomposition();
@@ -52,7 +54,8 @@ template<typename ModelType>
 std::shared_ptr<ModelType> performNondeterministicSparseBisimulationMinimization(std::shared_ptr<ModelType> model,
                                                                                  std::vector<std::shared_ptr<storm::logic::Formula const>> const& formulas,
                                                                                  storm::storage::BisimulationType type, bool graphPreserving = true,
-                                                                                 std::optional<double> const& tolerance = std::nullopt) {
+                                                                                 std::optional<double> const& tolerance = std::nullopt,
+                                                                                 bool allowMeasureDrivenInitialPartition = true) {
     using OptionsType = typename storm::storage::NondeterministicModelBisimulationDecomposition<ModelType>::Options;
     // Falls back to the general precision setting when the caller does not deliberately choose a tolerance;
     // may be reworked to require an explicit choice throughout the API in the future.
@@ -70,6 +73,7 @@ std::shared_ptr<ModelType> performNondeterministicSparseBisimulationMinimization
         options.setKeepRewards(true);
     }
     options.setType(type);
+    options.setAllowMeasureDrivenInitialPartition(allowMeasureDrivenInitialPartition);
 
     storm::storage::NondeterministicModelBisimulationDecomposition<ModelType> bisimulationDecomposition(*model, options);
     bisimulationDecomposition.computeBisimulationDecomposition();
@@ -80,7 +84,7 @@ template<typename ValueType>
 std::shared_ptr<storm::models::sparse::Model<ValueType>> performBisimulationMinimization(
     std::shared_ptr<storm::models::sparse::Model<ValueType>> const& model, std::vector<std::shared_ptr<storm::logic::Formula const>> const& formulas,
     storm::storage::BisimulationType type = storm::storage::BisimulationType::Strong, bool graphPreserving = true,
-    std::optional<double> const& tolerance = std::nullopt) {
+    std::optional<double> const& tolerance = std::nullopt, bool allowMeasureDrivenInitialPartition = true) {
     STORM_LOG_THROW(
         model->isOfType(storm::models::ModelType::Dtmc) || model->isOfType(storm::models::ModelType::Ctmc) || model->isOfType(storm::models::ModelType::Mdp),
         storm::exceptions::NotSupportedException, "Bisimulation minimization is currently only available for DTMCs, CTMCs and MDPs.");
@@ -90,13 +94,13 @@ std::shared_ptr<storm::models::sparse::Model<ValueType>> performBisimulationMini
 
     if (model->isOfType(storm::models::ModelType::Dtmc)) {
         return performDeterministicSparseBisimulationMinimization<storm::models::sparse::Dtmc<ValueType>>(
-            model->template as<storm::models::sparse::Dtmc<ValueType>>(), formulas, type, graphPreserving, tolerance);
+            model->template as<storm::models::sparse::Dtmc<ValueType>>(), formulas, type, graphPreserving, tolerance, allowMeasureDrivenInitialPartition);
     } else if (model->isOfType(storm::models::ModelType::Ctmc)) {
         return performDeterministicSparseBisimulationMinimization<storm::models::sparse::Ctmc<ValueType>>(
-            model->template as<storm::models::sparse::Ctmc<ValueType>>(), formulas, type, graphPreserving, tolerance);
+            model->template as<storm::models::sparse::Ctmc<ValueType>>(), formulas, type, graphPreserving, tolerance, allowMeasureDrivenInitialPartition);
     } else {
         return performNondeterministicSparseBisimulationMinimization<storm::models::sparse::Mdp<ValueType>>(
-            model->template as<storm::models::sparse::Mdp<ValueType>>(), formulas, type, graphPreserving, tolerance);
+            model->template as<storm::models::sparse::Mdp<ValueType>>(), formulas, type, graphPreserving, tolerance, allowMeasureDrivenInitialPartition);
     }
 }
 
diff --git a/src/storm/settings/modules/BisimulationSettings.cpp b/src/storm/settings/modules/BisimulationSettings.cpp
index 306ba2228..838d7ca48 100644
--- a/src/storm/settings/modules/BisimulationSettings.cpp
+++ b/src/storm/settings/modules/BisimulationSettings.cpp
@@ -22,6 +22,7 @@ const std::string BisimulationSettings::reuseOptionName = "reuse";
 const std::string BisimulationSettings::initialPartitionOptionName = "init";
 const std::string BisimulationSettings::refinementModeOptionName = "refine";
 const std::string BisimulationSettings::exactArithmeticDdOptionName = "ddexact";
+const std::string BisimulationSettings::noMeasureDrivenOptionName = "nomeasuredriven";
 
 BisimulationSettings::BisimulationSettings() : ModuleSettings(moduleName) {
     std::vector<std::string> types = {"strong", "weak"};
@@ -55,6 +56,10 @@ BisimulationSettings::BisimulationSettings() : ModuleSettings(moduleName) {
         storm::settings::OptionBuilder(moduleName, exactArithmeticDdOptionName, false, "Sets whether to use exact arithmetic in dd-based bisimulation.")
             .setIsAdvanced()
             .build());
+    this->addOption(storm::settings::OptionBuilder(moduleName, noMeasureDrivenOptionName, false,
+                                                   "Disables the measure-driven initial partition for sparse bisimulation minimization.")
+                        .setIsAdvanced()
+                        .build());
 
     std::vector<std::string> signatureModes = {"eager", "lazy"};
     this->addOption(storm::settings::OptionBuilder(moduleName, signatureModeOptionName, false, "Sets the signature computation mode.")
@@ -173,6 +178,10 @@ storm::dd::bisimulation::RefinementMode BisimulationSettings::getRefinementMode(
     return storm::dd::bisimulation::RefinementMode::Full;
 }
 
+bool BisimulationSettings::isMeasureDrivenPartitionDisabled() const {
+    return this->getOption(noMeasureDrivenOptionName).getHasOptionBeenSet();
+}
+
 bool BisimulationSettings::check() const {
     bool optionsSet = this->getOption(typeOptionName).getHasOptionBeenSet();
     STORM_LOG_WARN_COND(storm::settings::getModule<storm::settings::modules::GeneralSettings>().isBisimulationSet() || !optionsSet,
diff --git a/src/storm/settings/modules/BisimulationSettings.h b/src/storm/settings/modules/BisimulationSettings.h
index 435bef851..792aefc88 100644
--- a/src/storm/settings/modules/BisimulationSettings.h
+++ b/src/storm/settings/modules/BisimulationSettings.h
@@ -92,6 +92,12 @@ class BisimulationSettings : public ModuleSettings {
      */
     storm::dd::bisimulation::RefinementMode getRefinementMode() const;
 
+    /*!
+     * Retrieves whether the measure-driven initial partition is disabled for (sparse) bisimulation minimization.
+     * NOTE: only applies to sparse-based bisimulation.
+     */
+    bool isMeasureDrivenPartitionDisabled() const;
+
     virtual bool check() const override;
 
     // The name of the module.
@@ -109,6 +115,7 @@ class BisimulationSettings : public ModuleSettings {
     static const std::string refinementModeOptionName;
     static const std::string parallelismModeOptionName;
     static const std::string exactArithmeticDdOptionName;
+    static const std::string noMeasureDrivenOptionName;
 };
 }  // namespace modules
 }  // namespace settings
diff --git a/src/storm/storage/bisimulation/BisimulationDecomposition.cpp b/src/storm/storage/bisimulation/BisimulationDecomposition.cpp
index e82a60d4d..315559c78 100644
--- a/src/storm/storage/bisimulation/BisimulationDecomposition.cpp
+++ b/src/storm/storage/bisimulation/BisimulationDecomposition.cpp
@@ -128,21 +128,31 @@ void BisimulationDecomposition<ModelType, BlockDataType>::Options::checkAndSetMe
                 optimalityType = OptimizationDirection::Minimize;
             }
         }
+        formulaIsRewardObjective = true;
         newFormula = formula.asRewardOperatorFormula().getSubformula().asSharedPointer();
     }
 
     std::shared_ptr<storm::logic::Formula const> leftSubformula = std::make_shared<storm::logic::BooleanLiteralFormula>(true);
     std::shared_ptr<storm::logic::Formula const> rightSubformula;
     if (newFormula->isUntilFormula()) {
-        leftSubformula = newFormula->asUntilFormula().getLeftSubformula().asSharedPointer();
-        rightSubformula = newFormula->asUntilFormula().getRightSubformula().asSharedPointer();
-        if (leftSubformula->isInFragment(storm::logic::propositional()) && rightSubformula->isInFragment(storm::logic::propositional())) {
-            measureDrivenInitialPartition = true;
+        // UntilFormula is only ever probability-reachability; a reward objective can't reach this point.
+        if (!formulaIsRewardObjective) {
+            leftSubformula = newFormula->asUntilFormula().getLeftSubformula().asSharedPointer();
+            rightSubformula = newFormula->asUntilFormula().getRightSubformula().asSharedPointer();
+            if (leftSubformula->isInFragment(storm::logic::propositional()) && rightSubformula->isInFragment(storm::logic::propositional())) {
+                measureDrivenInitialPartition = allowMeasureDrivenInitialPartition;
+            }
         }
     } else if (newFormula->isEventuallyFormula()) {
-        rightSubformula = newFormula->asEventuallyFormula().getSubformula().asSharedPointer();
-        if (rightSubformula->isInFragment(storm::logic::propositional())) {
-            measureDrivenInitialPartition = true;
+        storm::logic::EventuallyFormula const& eventuallyFormula = newFormula->asEventuallyFormula();
+        rightSubformula = eventuallyFormula.getSubformula().asSharedPointer();
+
+        // Only the context matching the unwrapped operator (probability vs. reward) with default reward
+        // accumulation is supported by getStatesWithInfiniteReward/getStatesWithRewardZero.
+        bool const contextMatches =
+            formulaIsRewardObjective ? eventuallyFormula.isReachabilityRewardFormula() : eventuallyFormula.isReachabilityProbabilityFormula();
+        if (contextMatches && !eventuallyFormula.hasRewardAccumulation() && rightSubformula->isInFragment(storm::logic::propositional())) {
+            measureDrivenInitialPartition = allowMeasureDrivenInitialPartition;
         }
     }
 
@@ -354,17 +364,33 @@ void BisimulationDecomposition<ModelType, BlockDataType>::initializeLabelBasedPa
 }
 
 template<typename ModelType, typename BlockDataType>
-void BisimulationDecomposition<ModelType, BlockDataType>::initializeMeasureDrivenPartition() {
-    std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = this->getStatesWithProbability01();
+storm::storage::BitVector BisimulationDecomposition<ModelType, BlockDataType>::getStatesWithInfiniteReward() {
+    return this->getStatesWithProbability01().first;
+}
+
+template<typename ModelType, typename BlockDataType>
+storm::storage::BitVector BisimulationDecomposition<ModelType, BlockDataType>::getStatesWithRewardZero() {
+    return options.psiStates.value();
+}
 
+template<typename ModelType, typename BlockDataType>
+void BisimulationDecomposition<ModelType, BlockDataType>::initializeMeasureDrivenPartition() {
     std::optional<storm::storage::sparse::state_type> representativePsiState;
     if (!options.psiStates.value().empty()) {
         representativePsiState = *options.psiStates.value().begin();
     }
 
-    partition = storm::storage::bisimulation::Partition<BlockDataType>(
-        model.getNumberOfStates(), statesWithProbability01.first,
-        options.getBounded() || options.getKeepRewards() ? options.psiStates.value() : statesWithProbability01.second, representativePsiState);
+    if (options.formulaIsRewardObjective) {
+        // No meaningful "probability 1" counterpart for rewards; merge with the reward-zero states instead.
+        storm::storage::BitVector infinityStates = this->getStatesWithInfiniteReward();
+        storm::storage::BitVector rewardZeroStates = this->getStatesWithRewardZero();
+        partition = storm::storage::bisimulation::Partition<BlockDataType>(model.getNumberOfStates(), infinityStates, rewardZeroStates, representativePsiState);
+    } else {
+        std::pair<storm::storage::BitVector, storm::storage::BitVector> statesWithProbability01 = this->getStatesWithProbability01();
+        partition = storm::storage::bisimulation::Partition<BlockDataType>(model.getNumberOfStates(), statesWithProbability01.first,
+                                                                           options.getBounded() ? options.psiStates.value() : statesWithProbability01.second,
+                                                                           representativePsiState);
+    }
 
     // If the model has state rewards, we need to consider them, because otherwise reward properties are not
     // preserved.
diff --git a/src/storm/storage/bisimulation/BisimulationDecomposition.h b/src/storm/storage/bisimulation/BisimulationDecomposition.h
index 62f7d0b3f..35f8c459c 100644
--- a/src/storm/storage/bisimulation/BisimulationDecomposition.h
+++ b/src/storm/storage/bisimulation/BisimulationDecomposition.h
@@ -98,6 +98,15 @@ class BisimulationDecomposition : public Decomposition<StateBlock> {
             this->keepRewards = keepRewards;
         }
 
+        bool getAllowMeasureDrivenInitialPartition() const {
+            return this->allowMeasureDrivenInitialPartition;
+        }
+
+        // Disables the measure-driven initial partition even if it would otherwise be applicable.
+        void setAllowMeasureDrivenInitialPartition(bool value) {
+            this->allowMeasureDrivenInitialPartition = value;
+        }
+
         bool isOptimizationDirectionSet() const {
             return static_cast<bool>(optimalityType);
         }
@@ -122,6 +131,10 @@ class BisimulationDecomposition : public Decomposition<StateBlock> {
         std::optional<storm::storage::BitVector> phiStates;
         std::optional<storm::storage::BitVector> psiStates;
 
+        // Whether optimalityType refers to a reward objective (Rmin/Rmax) rather than a probability objective
+        // (Pmin/Pmax). Only meaningful together with measureDrivenInitialPartition.
+        bool formulaIsRewardObjective = false;
+
         /// An optional set of strings that indicate which of the atomic propositions of the model are to be
         /// respected and which may be ignored. If not given, all atomic propositions of the model are respected.
         std::optional<std::set<std::string>> respectedAtomicPropositions;
@@ -143,6 +156,9 @@ class BisimulationDecomposition : public Decomposition<StateBlock> {
         /// when computing strong bisimulation equivalence.
         bool bounded = false;
 
+        /// Whether a measure-driven initial partition may be used (if applicable).
+        bool allowMeasureDrivenInitialPartition = true;
+
         /// A flag that indicates whether discounted properties are to be preserved. This may only be set to true
         /// when computing strong bisimulation equivalence.
         bool discounted = false;
@@ -262,6 +278,21 @@ class BisimulationDecomposition : public Decomposition<StateBlock> {
      */
     virtual std::pair<storm::storage::BitVector, storm::storage::BitVector> getStatesWithProbability01() = 0;
 
+    /*!
+     * Computes the states with infinite expected reward until psi. Used instead of getStatesWithProbability01()
+     * for the measure-driven initial partition of a reward objective. Default implementation delegates to
+     * getStatesWithProbability01(); nondeterministic models override this (see
+     * NondeterministicModelBisimulationDecomposition::getStatesWithInfiniteReward).
+     */
+    virtual storm::storage::BitVector getStatesWithInfiniteReward();
+
+    /*!
+     * Computes the states known a priori to have expected reward 0 until psi (a superset of the psi states
+     * themselves). Default implementation just returns the psi states; nondeterministic models can compute a
+     * larger set (see NondeterministicModelBisimulationDecomposition::getStatesWithRewardZero).
+     */
+    virtual storm::storage::BitVector getStatesWithRewardZero();
+
     /*!
      * Splits the initial partition based on the (unique) reward model of the current model.
      */
diff --git a/src/storm/storage/bisimulation/DeterministicBlockData.cpp b/src/storm/storage/bisimulation/DeterministicBlockData.cpp
index b7e2a8e2b..9530d6aec 100644
--- a/src/storm/storage/bisimulation/DeterministicBlockData.cpp
+++ b/src/storm/storage/bisimulation/DeterministicBlockData.cpp
@@ -93,6 +93,22 @@ storm::storage::sparse::state_type DeterministicBlockData::representativeState()
     return valRepresentativeState.value();
 }
 
+void DeterministicBlockData::setProb0(bool value) {
+    setFlag(PROB0_FLAG, value);
+}
+
+bool DeterministicBlockData::prob0() const {
+    return getFlag(PROB0_FLAG);
+}
+
+void DeterministicBlockData::setProb1(bool value) {
+    setFlag(PROB1_FLAG, value);
+}
+
+bool DeterministicBlockData::prob1() const {
+    return getFlag(PROB1_FLAG);
+}
+
 bool DeterministicBlockData::needsRefinement() const {
     return getFlag(REFINEMENT_FLAG);
 }
diff --git a/src/storm/storage/bisimulation/DeterministicBlockData.h b/src/storm/storage/bisimulation/DeterministicBlockData.h
index 84f274436..860da9879 100644
--- a/src/storm/storage/bisimulation/DeterministicBlockData.h
+++ b/src/storm/storage/bisimulation/DeterministicBlockData.h
@@ -65,6 +65,18 @@ class DeterministicBlockData {
     // Retrieves the representative state for this block.
     storm::storage::sparse::state_type representativeState() const;
 
+    // Sets whether this is the prob0 block of a measure-driven initial partition.
+    void setProb0(bool value = true);
+
+    // Retrieves whether this is the prob0 block of a measure-driven initial partition.
+    bool prob0() const;
+
+    // Sets whether this is the prob1 block of a measure-driven initial partition.
+    void setProb1(bool value = true);
+
+    // Retrieves whether this is the prob1 block of a measure-driven initial partition.
+    bool prob1() const;
+
     friend std::ostream& operator<<(std::ostream& out, DeterministicBlockData const& data);
 
    public:
@@ -82,6 +94,8 @@ class DeterministicBlockData {
     static constexpr uint64_t REFINEMENT_FLAG = 1ull << 1;
     static constexpr uint64_t ABSORBING_FLAG = 1ull << 2;
     static constexpr uint64_t REWARD_FLAG = 1ull << 3;
+    static constexpr uint64_t PROB0_FLAG = 1ull << 4;
+    static constexpr uint64_t PROB1_FLAG = 1ull << 5;
     uint8_t flags;
 
     // An optional representative state for the block. If this is set, this state is used to derive the
diff --git a/src/storm/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp b/src/storm/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp
index b422d19c1..30765ab95 100644
--- a/src/storm/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp
+++ b/src/storm/storage/bisimulation/DeterministicModelBisimulationDecomposition.cpp
@@ -601,6 +601,16 @@ void DeterministicModelBisimulationDecomposition<ModelType>::buildQuotient() {
         newLabeling.addLabel(ap);
     }
 
+    // Mark synthetic blocks produced by a measure-driven initial partition, so they're distinguishable from
+    // blocks that correspond to a single original state.
+    bool const addMeasureDrivenLabels = this->options.measureDrivenInitialPartition;
+    std::string const firstBlockLabel = this->options.formulaIsRewardObjective ? "__rewinf__" : "__prob0__";
+    std::string const secondBlockLabel = this->options.formulaIsRewardObjective ? "__rew0__" : "__prob1__";
+    if (addMeasureDrivenLabels) {
+        newLabeling.addLabel(firstBlockLabel);
+        newLabeling.addLabel(secondBlockLabel);
+    }
+
     // If the model had state rewards, we need to build the state rewards for the quotient as well.
     std::optional<std::vector<ValueType>> stateRewards;
     if (this->options.getKeepRewards() && this->model.hasRewardModel()) {
@@ -644,6 +654,14 @@ void DeterministicModelBisimulationDecomposition<ModelType>::buildQuotient() {
                     newLabeling.addLabelToState(ap, blockIndex);
                 }
             }
+
+            if (addMeasureDrivenLabels) {
+                if (oldBlock.data().prob0()) {
+                    newLabeling.addLabelToState(firstBlockLabel, blockIndex);
+                } else if (oldBlock.data().prob1()) {
+                    newLabeling.addLabelToState(secondBlockLabel, blockIndex);
+                }
+            }
         } else {
             // Compute the outgoing transitions of the block.
             std::map<storm::storage::sparse::state_type, ValueType> blockProbability;
diff --git a/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp b/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp
index e70eabdda..a76e34bcb 100644
--- a/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp
+++ b/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.cpp
@@ -37,6 +37,45 @@ std::pair<storm::storage::BitVector, storm::storage::BitVector> Nondeterministic
     }
 }
 
+template<typename ModelType>
+storm::storage::BitVector NondeterministicModelBisimulationDecomposition<ModelType>::getStatesWithInfiniteReward() {
+    STORM_LOG_THROW(this->options.isOptimizationDirectionSet(), storm::exceptions::IllegalFunctionCallException,
+                    "Can only compute states with infinite reward with an optimization direction (min/max).");
+    // Rmin is infinite iff no scheduler reaches psi a.s. (not in the existential Prob1E set); Rmax is infinite
+    // iff not every scheduler does (not in the universal Prob1A set) - the opposite of the reward's own
+    // direction. Cf. SparseMdpPrctlHelper::computeQualitativeStateSetsReachabilityRewards.
+    if (this->options.getOptimizationDirection() == OptimizationDirection::Minimize) {
+        return ~storm::utility::graph::performProb1E(this->model.getTransitionMatrix(), this->model.getTransitionMatrix().getRowGroupIndices(),
+                                                     this->model.getBackwardTransitions(), this->options.phiStates.value(), this->options.psiStates.value());
+    } else {
+        return ~storm::utility::graph::performProb1A(this->model.getTransitionMatrix(), this->model.getTransitionMatrix().getRowGroupIndices(),
+                                                     this->model.getBackwardTransitions(), this->options.phiStates.value(), this->options.psiStates.value());
+    }
+}
+
+template<typename ModelType>
+storm::storage::BitVector NondeterministicModelBisimulationDecomposition<ModelType>::getStatesWithRewardZero() {
+    if (!this->model.hasRewardModel()) {
+        return this->options.psiStates.value();
+    }
+    STORM_LOG_THROW(this->options.isOptimizationDirectionSet(), storm::exceptions::IllegalFunctionCallException,
+                    "Can only compute states with reward zero with an optimization direction (min/max).");
+    auto const& rewardModel = this->model.getUniqueRewardModel();
+    storm::storage::BitVector trueStates(this->model.getNumberOfStates(), true);
+
+    // For Rmin: reward 0 if some scheduler reaches psi using only zero-reward choices (existential). For Rmax:
+    // reward 0 only if every scheduler staying within zero-reward states reaches psi (universal).
+    if (this->options.getOptimizationDirection() == OptimizationDirection::Minimize) {
+        storm::storage::BitVector zeroRewardChoices = rewardModel.getChoicesWithZeroReward(this->model.getTransitionMatrix());
+        return storm::utility::graph::performProb1E(this->model.getTransitionMatrix(), this->model.getTransitionMatrix().getRowGroupIndices(),
+                                                    this->model.getBackwardTransitions(), trueStates, this->options.psiStates.value(), zeroRewardChoices);
+    } else {
+        storm::storage::BitVector zeroRewardStates = rewardModel.getStatesWithZeroReward(this->model.getTransitionMatrix());
+        return storm::utility::graph::performProb1A(this->model.getTransitionMatrix(), this->model.getTransitionMatrix().getRowGroupIndices(),
+                                                    this->model.getBackwardTransitions(), zeroRewardStates, this->options.psiStates.value());
+    }
+}
+
 template<typename ModelType>
 void NondeterministicModelBisimulationDecomposition<ModelType>::initialize() {
     this->createChoiceToStateMapping();
@@ -121,6 +160,16 @@ void NondeterministicModelBisimulationDecomposition<ModelType>::buildQuotient()
         newLabeling.addLabel(ap);
     }
 
+    // Mark synthetic blocks produced by a measure-driven initial partition, so they're distinguishable from
+    // blocks that correspond to a single original state.
+    bool const addMeasureDrivenLabels = this->options.measureDrivenInitialPartition;
+    std::string const firstBlockLabel = this->options.formulaIsRewardObjective ? "__rewinf__" : "__prob0__";
+    std::string const secondBlockLabel = this->options.formulaIsRewardObjective ? "__rew0__" : "__prob1__";
+    if (addMeasureDrivenLabels) {
+        newLabeling.addLabel(firstBlockLabel);
+        newLabeling.addLabel(secondBlockLabel);
+    }
+
     // If the model had state (action) rewards, we need to build the state rewards for the quotient as well.
     std::optional<std::vector<ValueType>> stateRewards;
     std::optional<std::vector<ValueType>> stateActionRewards;
@@ -169,6 +218,14 @@ void NondeterministicModelBisimulationDecomposition<ModelType>::buildQuotient()
                     newLabeling.addLabelToState(ap, blockIndex);
                 }
             }
+
+            if (addMeasureDrivenLabels) {
+                if (oldBlock.data().prob0()) {
+                    newLabeling.addLabelToState(firstBlockLabel, blockIndex);
+                } else if (oldBlock.data().prob1()) {
+                    newLabeling.addLabelToState(secondBlockLabel, blockIndex);
+                }
+            }
         } else {
             // Add the outgoing choices of the block.
             for (uint_fast64_t choice = nondeterministicChoiceIndices[representativeState]; choice < nondeterministicChoiceIndices[representativeState + 1];
diff --git a/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h b/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h
index cffaf8da2..8f1e28205 100644
--- a/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h
+++ b/src/storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h
@@ -30,6 +30,10 @@ class NondeterministicModelBisimulationDecomposition : public BisimulationDecomp
    protected:
     virtual std::pair<storm::storage::BitVector, storm::storage::BitVector> getStatesWithProbability01() override;
 
+    virtual storm::storage::BitVector getStatesWithInfiniteReward() override;
+
+    virtual storm::storage::BitVector getStatesWithRewardZero() override;
+
     virtual void buildQuotient() override;
 
     virtual void refinePartitionBasedOnSplitter(bisimulation::Block<BlockDataType>& splitter,
diff --git a/src/storm/storage/bisimulation/Partition.cpp b/src/storm/storage/bisimulation/Partition.cpp
index 798168562..08324cbb6 100644
--- a/src/storm/storage/bisimulation/Partition.cpp
+++ b/src/storm/storage/bisimulation/Partition.cpp
@@ -40,6 +40,7 @@ Partition<DataType>::Partition(std::size_t numberOfStates, storm::storage::BitVe
             ++position;
         }
         firstBlock->data().setAbsorbing(true);
+        firstBlock->data().setProb0(true);
     }
 
     if (!prob1States.empty()) {
@@ -53,6 +54,7 @@ Partition<DataType>::Partition(std::size_t numberOfStates, storm::storage::BitVe
             ++position;
         }
         secondBlock->data().setAbsorbing(true);
+        secondBlock->data().setProb1(true);
         secondBlock->data().setRepresentativeState(representativeProb1State.value());
     }
 
diff --git a/src/test/storm/storage/NondeterministicModelBisimulationDecompositionTest.cpp b/src/test/storm/storage/NondeterministicModelBisimulationDecompositionTest.cpp
index 48c8cf888..50402a1a2 100644
--- a/src/test/storm/storage/NondeterministicModelBisimulationDecompositionTest.cpp
+++ b/src/test/storm/storage/NondeterministicModelBisimulationDecompositionTest.cpp
@@ -3,15 +3,54 @@
 
 #include "storm-parsers/parser/FormulaParser.h"
 #include "storm-parsers/parser/PrismParser.h"
+#include "storm/api/storm.h"
 #include "storm/builder/ExplicitModelBuilder.h"
+#include "storm/modelchecker/results/ExplicitQuantitativeCheckResult.h"
 #include "storm/models/sparse/Mdp.h"
 #include "storm/models/sparse/StandardRewardModel.h"
 #include "storm/storage/bisimulation/NondeterministicModelBisimulationDecomposition.h"
 
 namespace {
 static constexpr double DefaultTestTolerance = 1e-6;
+
+// Model checks propertyString on both the original model and its bisimulation quotient; returns (ground truth,
+// quotient value) for the initial state.
+std::pair<double, double> computeGroundTruthAndBisimulationResult(std::string const& modelFile, std::string const& propertyString) {
+    storm::prism::Program program = storm::parser::PrismParser::parse(modelFile);
+
+    storm::parser::FormulaParser formulaParser;
+    std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString(propertyString);
+    std::vector<std::shared_ptr<storm::logic::Formula const>> formulas = {formula};
+
+    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::api::buildSparseModel<double>(program, formulas);
+    std::shared_ptr<storm::models::sparse::Mdp<double>> mdp = model->as<storm::models::sparse::Mdp<double>>();
+
+    double groundTruth = storm::api::verifyWithSparseEngine(model, storm::api::createTask<double>(formula, true))
+                             ->template asExplicitQuantitativeCheckResult<double>()[*model->getInitialStates().begin()];
+
+    typename storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>>::Options options(*mdp, *formula,
+                                                                                                                                 DefaultTestTolerance);
+    storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>> bisim(*mdp, options);
+    bisim.computeBisimulationDecomposition();
+    std::shared_ptr<storm::models::sparse::Model<double>> quotient = bisim.getQuotient();
+
+    double quotientValue = storm::api::verifyWithSparseEngine(quotient, storm::api::createTask<double>(formula, true))
+                               ->template asExplicitQuantitativeCheckResult<double>()[*quotient->getInitialStates().begin()];
+
+    return {groundTruth, quotientValue};
 }
 
+// EXPECT_NEAR fails on infinities (inf - inf = NaN). Mirrors expectVectorNear in storm-dft/bdd/TestBdd.cpp.
+void expectRewardNear(double expected, double actual, double precision = 1e-6) {
+    if (std::isinf(expected)) {
+        EXPECT_EQ(expected, actual);
+    } else {
+        EXPECT_NEAR(expected, actual, precision);
+    }
+}
+
+}  // namespace
+
 TEST(NondeterministicModelBisimulationDecomposition, TwoDice) {
 #ifndef STORM_HAVE_Z3
     GTEST_SKIP() << "Z3 not available.";
@@ -65,3 +104,85 @@ TEST(NondeterministicModelBisimulationDecomposition, TwoDice) {
     EXPECT_EQ(26ul, result->getNumberOfTransitions());
     EXPECT_EQ(14ul, result->as<storm::models::sparse::Mdp<double>>()->getNumberOfChoices());
 }
+
+// Regression test: measure-driven partition used to merge all probability-1 states into one block regardless of
+// their (differing, finite) expected reward, turning a finite Rmin into "infinity" in the quotient.
+TEST(NondeterministicModelBisimulationDecomposition, MeasureDrivenRewardUnsound) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    storm::prism::Program program = storm::parser::PrismParser::parse(STORM_TEST_RESOURCES_DIR "/mdp/tiny_rewards2.nm");
+
+    storm::parser::FormulaParser formulaParser;
+    std::shared_ptr<storm::logic::Formula const> formula = formulaParser.parseSingleFormulaFromString("Rmin=? [F \"goal\"]");
+    std::vector<std::shared_ptr<storm::logic::Formula const>> formulas = {formula};
+
+    std::shared_ptr<storm::models::sparse::Model<double>> model = storm::api::buildSparseModel<double>(program, formulas);
+    ASSERT_EQ(model->getType(), storm::models::ModelType::Mdp);
+    std::shared_ptr<storm::models::sparse::Mdp<double>> mdp = model->as<storm::models::sparse::Mdp<double>>();
+    size_t initialState = *model->getInitialStates().begin();
+
+    // Ground truth: model check the original, unreduced model.
+    std::unique_ptr<storm::modelchecker::CheckResult> groundTruthCheckResult =
+        storm::api::verifyWithSparseEngine(model, storm::api::createTask<double>(formula, true));
+    double groundTruthValue = groundTruthCheckResult->asExplicitQuantitativeCheckResult<double>()[initialState];
+    // Sanity check: the correct result is finite (a scheduler exists that reaches "goal" almost surely).
+    ASSERT_LT(groundTruthValue, storm::utility::infinity<double>());
+
+    typename storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>>::Options options(*mdp, *formula,
+                                                                                                                                 DefaultTestTolerance);
+    storm::storage::NondeterministicModelBisimulationDecomposition<storm::models::sparse::Mdp<double>> bisim(*mdp, options);
+    ASSERT_NO_THROW(bisim.computeBisimulationDecomposition());
+    std::shared_ptr<storm::models::sparse::Model<double>> quotient;
+    ASSERT_NO_THROW(quotient = bisim.getQuotient());
+    size_t quotientInitialState = *quotient->getInitialStates().begin();
+
+    std::unique_ptr<storm::modelchecker::CheckResult> quotientCheckResult =
+        storm::api::verifyWithSparseEngine(quotient, storm::api::createTask<double>(formula, true));
+    double quotientValue = quotientCheckResult->asExplicitQuantitativeCheckResult<double>()[quotientInitialState];
+
+    // The quotient must preserve the value of the formula it was built for.
+    EXPECT_NEAR(groundTruthValue, quotientValue, 1e-6);
+}
+
+// Model from https://github.com/moves-rwth/storm/issues/683: action a3 (at s=2) loops back towards "goal" (finite
+// Rmin), while a4 self-loops at s=2 forever (infinite Rmax). Exercises the Prob1E/Prob1A asymmetry.
+TEST(NondeterministicModelBisimulationDecomposition, TinyRewards3Issue683) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    auto [groundTruthMin, quotientMin] = computeGroundTruthAndBisimulationResult(STORM_TEST_RESOURCES_DIR "/mdp/tiny_rewards3.nm", "Rmin=? [F \"goal\"]");
+    ASSERT_LT(groundTruthMin, storm::utility::infinity<double>());
+    expectRewardNear(groundTruthMin, quotientMin);
+
+    auto [groundTruthMax, quotientMax] = computeGroundTruthAndBisimulationResult(STORM_TEST_RESOURCES_DIR "/mdp/tiny_rewards3.nm", "Rmax=? [F \"goal\"]");
+    ASSERT_TRUE(std::isinf(groundTruthMax));
+    expectRewardNear(groundTruthMax, quotientMax);
+}
+
+// s=0 can self-loop forever ("s=0 -> true"), avoiding "target" without incurring reward. Rmax is still infinite
+// by convention, even though the self-loop's own reward is 0.
+TEST(NondeterministicModelBisimulationDecomposition, TinyRewardsInfinity) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    auto [groundTruthMin, quotientMin] = computeGroundTruthAndBisimulationResult(STORM_TEST_RESOURCES_DIR "/mdp/tiny_rewards.nm", "Rmin=? [F \"target\"]");
+    ASSERT_LT(groundTruthMin, storm::utility::infinity<double>());
+    expectRewardNear(groundTruthMin, quotientMin);
+
+    auto [groundTruthMax, quotientMax] = computeGroundTruthAndBisimulationResult(STORM_TEST_RESOURCES_DIR "/mdp/tiny_rewards.nm", "Rmax=? [F \"target\"]");
+    ASSERT_TRUE(std::isinf(groundTruthMax));
+    expectRewardNear(groundTruthMax, quotientMax);
+}
+
+// Coverage on a larger, non-artificial model: "coinflips" reward on two_dice.nm (Rmin == Rmax, no nondeterminism).
+TEST(NondeterministicModelBisimulationDecomposition, TwoDiceExpectedCoinFlips) {
+#ifndef STORM_HAVE_Z3
+    GTEST_SKIP() << "Z3 not available.";
+#endif
+    auto [groundTruthMin, quotientMin] = computeGroundTruthAndBisimulationResult(STORM_TEST_RESOURCES_DIR "/mdp/two_dice.nm", "Rmin=? [F \"done\"]");
+    expectRewardNear(groundTruthMin, quotientMin);
+
+    auto [groundTruthMax, quotientMax] = computeGroundTruthAndBisimulationResult(STORM_TEST_RESOURCES_DIR "/mdp/two_dice.nm", "Rmax=? [F \"done\"]");
+    expectRewardNear(groundTruthMax, quotientMax);
+}