diff --git a/stan/math/rev/functor/cvodes_integrator_adjoint.hpp b/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
index d8e030355ee..24ae975debc 100644
--- a/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
+++ b/stan/math/rev/functor/cvodes_integrator_adjoint.hpp
@@ -641,6 +641,9 @@ class cvodes_integrator_adjoint_vari : public vari_base {
                      f_y_t_vars.size(), "states", N_);
     f_y_t_vars.adj() = -Eigen::Map<Eigen::VectorXd>(NV_DATA_S(yB), N_);
     grad();
+    // CVODES reuses this output buffer; overwrite its previous contents before
+    // accumulating the current quadrature derivatives.
+    N_VConst(0.0, qBdot);
     math::apply(
         [&qBdot](auto&&... args) {
           accumulate_adjoints(NV_DATA_S(qBdot), args...);
diff --git a/test/unit/math/rev/functor/cos_ode_typed_test.cpp b/test/unit/math/rev/functor/cos_ode_typed_test.cpp
index 327842289ee..f945120312d 100644
--- a/test/unit/math/rev/functor/cos_ode_typed_test.cpp
+++ b/test/unit/math/rev/functor/cos_ode_typed_test.cpp
@@ -5,6 +5,7 @@
 #include <test/unit/math/rev/functor/test_fixture_ode.hpp>
 #include <test/unit/math/rev/functor/test_fixture_ode_cos_scalar.hpp>
 #include <test/unit/math/rev/functor/ode_test_functors.hpp>
+#include <cmath>
 
 /**
  *
@@ -96,3 +97,32 @@ REGISTER_TYPED_TEST_SUITE_P(cos_arg_test, y0_error, t0_error, ts_error,
                             one_arg_error, two_arg_error, rhs_wrong_size_error,
                             error_name, tol_error, value, grad, tol_grad);
 INSTANTIATE_TYPED_TEST_SUITE_P(StanOde, cos_arg_test, cos_arg_test_types);
+
+// CVODES buffer reuse: SUNDIALS 6.1.1, src/cvodes/cvodes.c:6522-6536.
+TEST_F(AgradRev, adjoint_quadrature_buffer_reuse) {
+  const std::vector<double> times{0.45, 1.1};
+  constexpr double rate = 1.5;
+  for (int iteration = 0; iteration < 5; ++iteration) {
+    SCOPED_TRACE(iteration);
+    stan::math::nested_rev_autodiff nested;
+    const Eigen::VectorXd initial = Eigen::VectorXd::Zero(1);
+    stan::math::var scalar = 0.75;
+    std::vector<stan::math::var> parameters{0.25, 0.5};
+    ode_adjoint_functor solver;
+    auto output = solver(stan::test::Cos2Arg(), initial, 0.0, times, nullptr,
+                         scalar, parameters);
+    for (int index : {0, 1, 0, 1}) {
+      SCOPED_TRACE(index);
+      nested.set_zero_all_adjoints();
+      ASSERT_NO_THROW(output[index](0).grad());
+      const double phase = rate * times[index];
+      const double expected
+          = (phase * std::cos(phase) - std::sin(phase)) / (rate * rate);
+      EXPECT_NEAR(output[index](0).val(), std::sin(phase) / rate, 1e-8);
+      EXPECT_NEAR(scalar.adj(), expected, 1e-8);
+      for (const auto& parameter : parameters) {
+        EXPECT_NEAR(parameter.adj(), expected, 1e-8);
+      }
+    }
+  }
+}