diff --git a/src/api/api_qe.cpp b/src/api/api_qe.cpp
index 94067135f..5424c0330 100644
--- a/src/api/api_qe.cpp
+++ b/src/api/api_qe.cpp
@@ -26,6 +26,7 @@ Notes:
 #include "api/api_ast_map.h"
 #include "api/api_ast_vector.h"
 #include "qe/lite/qe_lite_tactic.h"
+#include "qe/qe_mbp.h"
 #include "muz/spacer/spacer_util.h"
 
 extern "C"
@@ -104,6 +105,46 @@ extern "C"
         Z3_CATCH_RETURN(nullptr);
     }
 
+    Z3_ast Z3_API Z3_qe_model_project_with_witness (Z3_context c,
+                                              Z3_model mdl,
+                                              unsigned num_bounds,
+                                              Z3_app const bound[],
+                                              Z3_ast body,
+                                              Z3_ast_map map)
+    {
+        Z3_TRY;
+        LOG_Z3_qe_model_project_with_witness (c, mdl, num_bounds, bound, body, map);
+        RESET_ERROR_CODE();
+
+        ast_manager& m = mk_c(c)->m();
+        app_ref_vector vars(m);
+        if (!to_apps(num_bounds, bound, vars)) {
+            RETURN_Z3(nullptr);
+        }
+
+        expr_ref_vector fmls(m);
+        fmls.push_back(to_expr(body));
+        model_ref model (to_model_ref (mdl));
+        vector<mbp::def> defs;
+        qe::mbproj proj(m);
+
+        proj(true, vars, *model, fmls, &defs);
+        expr_ref result (m);
+        result = m.mk_and(fmls);
+        mk_c(c)->save_ast_trail(result);
+
+        obj_map<ast, ast*> &map_z3 = to_ast_map_ref(map);
+
+        for (auto& [v, t] : defs) {
+            m.inc_ref(v);
+            m.inc_ref(t);
+            map_z3.insert(v, t);            
+        }
+
+        return of_expr (result);
+        Z3_CATCH_RETURN(nullptr);
+    }
+
     Z3_ast Z3_API Z3_model_extrapolate (Z3_context c,
                                         Z3_model m,
                                         Z3_ast fml)
diff --git a/src/api/python/z3/z3.py b/src/api/python/z3/z3.py
index 7aac297a1..bd59c9917 100644
--- a/src/api/python/z3/z3.py
+++ b/src/api/python/z3/z3.py
@@ -6732,6 +6732,30 @@ class ModelRef(Z3PPObject):
         model = Z3_model_translate(self.ctx.ref(), self.model, target.ref())
         return ModelRef(model, target)
 
+    def project(self, vars, fml):
+        """Perform model-based projection on fml with respect to vars.
+        Assume that the model satisfies fml. Then compute a projection fml_p, such
+        that vars do not occur free in fml_p, fml_p is true in the model and
+        fml_p => exists vars . fml
+        """
+        ctx = self.ctx.ref()
+        _vars = (Ast * len(vars))()
+        for i in range(len(vars)):
+            _vars[i] = vars[i].as_ast()
+        return _to_expr_ref(Z3_qe_model_project(ctx, self.model, len(vars), _vars, fml.ast), self.ctx)
+
+    def project_with_witness(self, vars, fml):
+        """Perform model-based projection, but also include realizer terms for the projected variables"""
+        ctx = self.ctx.ref()
+        _vars = (Ast * len(vars))()
+        for i in range(len(vars)):
+            _vars[i] = vars[i].as_ast()
+        defs = AstMap()
+        result = Z3_qe_model_project_with_witness(ctx, self.model, len(vars), _vars, fml.ast, defs.map)
+        result = _to_expr_ref(result, self.ctx)
+        return result, defs
+
+
     def __copy__(self):
         return self.translate(self.ctx)
 
@@ -6739,9 +6763,12 @@ class ModelRef(Z3PPObject):
         return self.translate(self.ctx)
 
 
-def Model(ctx=None):
+def Model(ctx=None, eval = {}):
     ctx = _get_ctx(ctx)
-    return ModelRef(Z3_mk_model(ctx.ref()), ctx)
+    mdl = ModelRef(Z3_mk_model(ctx.ref()), ctx)
+    for k, v in eval.items():
+        mdl.update_value(k, v)
+    return mdl
 
 
 def is_as_array(n):
diff --git a/src/api/z3_spacer.h b/src/api/z3_spacer.h
index 1f7a7ef34..81d419275 100644
--- a/src/api/z3_spacer.h
+++ b/src/api/z3_spacer.h
@@ -112,6 +112,22 @@ extern "C" {
        Z3_ast body,
        Z3_ast_map map);
 
+    /**
+       \brief Project with witness extraction.
+
+       The returned map contains a binding of variables to terms that such that when the binding
+       is used for the formula, it remains true within the model.
+
+       def_API('Z3_qe_model_project_with_witness', AST, (_in(CONTEXT), _in(MODEL), _in(UINT), _in_array(2, APP), _in(AST), _in(AST_MAP)))
+    */
+    Z3_ast Z3_API Z3_qe_model_project_with_witness
+      (Z3_context c,
+       Z3_model m,
+       unsigned num_bounds,
+       Z3_app const bound[],
+       Z3_ast body,
+       Z3_ast_map map);
+
     /**
        \brief Extrapolates a model of a formula
 
diff --git a/src/qe/mbp/mbp_plugin.h b/src/qe/mbp/mbp_plugin.h
index c270489e8..5cf91098b 100644
--- a/src/qe/mbp/mbp_plugin.h
+++ b/src/qe/mbp/mbp_plugin.h
@@ -28,7 +28,7 @@ Revision History:
 
 namespace mbp {
 
-    struct cant_project {};
+    struct cant_project : public std::exception {};
 
     struct def {
         expr_ref var, term;
diff --git a/src/qe/qe_mbp.cpp b/src/qe/qe_mbp.cpp
index 5f5752187..27751adb9 100644
--- a/src/qe/qe_mbp.cpp
+++ b/src/qe/qe_mbp.cpp
@@ -418,12 +418,12 @@ public:
         e = mk_and(fmls);
         return any_of(subterms::all(e), [&](expr* c) { return seq.is_char(c) || seq.is_seq(c); });
     }
-    void operator()(bool force_elim, app_ref_vector& vars, model& model, expr_ref_vector& fmls) {
+    void operator()(bool force_elim, app_ref_vector& vars, model& model, expr_ref_vector& fmls, vector<mbp::def>* defs = nullptr) {
             //don't use mbp_qel on some theories where model evaluation is
             //incomplete This is not a limitation of qel. Fix this either by
             //making mbp choices less dependent on the model evaluation methods
             //or fix theory rewriters to make terms evaluation complete
-            if (m_use_qel && !has_unsupported_th(fmls)) {
+            if (m_use_qel && !has_unsupported_th(fmls) && !defs) {
                 bool dsub = m_dont_sub;
                 m_dont_sub = !force_elim;
                 expr_ref fml(m);
@@ -434,11 +434,11 @@ public:
                 m_dont_sub = dsub;
             }
             else {
-                mbp(force_elim, vars, model, fmls);
+                mbp(force_elim, vars, model, fmls, defs);
             }
         }
 
-    void mbp(bool force_elim, app_ref_vector& vars, model& model, expr_ref_vector& fmls) {
+    void mbp(bool force_elim, app_ref_vector& vars, model& model, expr_ref_vector& fmls, vector<mbp::def>* defs) {
         SASSERT(validate_model(model, fmls));
         expr_ref val(m), tmp(m);
         app_ref var(m);
@@ -451,10 +451,15 @@ public:
             app_ref_vector new_vars(m);
             progress = false;
             for (mbp::project_plugin* p : m_plugins) {
-                if (p)
+                if (defs && p) {
+                    unsigned sz = defs->size();
+                    p->project(model, vars, fmls, *defs);
+                    progress |= sz < defs->size();
+                }
+                else if (p)
                     (*p)(model, vars, fmls);
             }
-            while (!vars.empty() && !fmls.empty() && m.limit().inc()) {
+            while (!vars.empty() && !fmls.empty() && !defs && m.limit().inc()) {
                 var = vars.back();
                 vars.pop_back();
                 mbp::project_plugin* p = get_plugin(var);
@@ -471,6 +476,8 @@ public:
                 expr_safe_replace sub(m);
                 val = model(var);
                 sub.insert(var, val);
+                if (defs)
+                    defs->push_back(mbp::def(expr_ref(var, m), val));
                 unsigned j = 0;
                 for (expr* f : fmls) {
                     sub(f, tmp);
@@ -562,7 +569,7 @@ public:
             expr_ref_vector fmls(m);
             flatten_and(fml, fmls);
 
-            mbp(false, other_vars, mdl, fmls);
+            mbp(false, other_vars, mdl, fmls, nullptr);
             fml = mk_and(fmls);
             m_rw(fml);
 
@@ -704,9 +711,9 @@ void mbproj::get_param_descrs(param_descrs& r) {
     r.insert("use_qel", CPK_BOOL, "(default: true) use egraph based QEL");
 }
 
-void mbproj::operator()(bool force_elim, app_ref_vector& vars, model& mdl, expr_ref_vector& fmls) {
+void mbproj::operator()(bool force_elim, app_ref_vector& vars, model& mdl, expr_ref_vector& fmls, vector<mbp::def>* defs) {
     scoped_no_proof _sp(fmls.get_manager());
-    (*m_impl)(force_elim, vars, mdl, fmls);
+    (*m_impl)(force_elim, vars, mdl, fmls, defs);
 }
 
 void mbproj::spacer(app_ref_vector& vars, model& mdl, expr_ref& fml) {
diff --git a/src/qe/qe_mbp.h b/src/qe/qe_mbp.h
index daa65de7a..3d85d3e88 100644
--- a/src/qe/qe_mbp.h
+++ b/src/qe/qe_mbp.h
@@ -24,6 +24,7 @@ Revision History:
 #include "util/params.h"
 #include "model/model.h"
 #include "math/simplex/model_based_opt.h"
+#include "qe/mbp/mbp_plugin.h"
 
 
 namespace qe {
@@ -45,7 +46,7 @@ namespace qe {
            Apply model-based qe on constants provided as vector of variables.
            Return the updated formula and updated set of variables that were not eliminated.
         */
-        void operator()(bool force_elim, app_ref_vector& vars, model& mdl, expr_ref_vector& fmls);
+        void operator()(bool force_elim, app_ref_vector& vars, model& mdl, expr_ref_vector& fmls, vector<mbp::def>* defs = nullptr);
 
         /**
            \brief
@@ -70,3 +71,4 @@ namespace qe {
     };
 }
 
+