diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst
index 9e5f0d819..204d262e3 100644
--- a/Docs/sphinx_doc/theory/Buoyancy.rst
+++ b/Docs/sphinx_doc/theory/Buoyancy.rst
@@ -7,6 +7,29 @@
 
 .. _Buoyancy:
 
+Density of the mixture
+========================
+
+The total density in a given cell is given by
+
+.. math::
+    \rho = \frac{m}{V} = \frac{m_a + m_v + m_c + m_p}{V},
+
+where :math:`m_a` is the mass of dry air, :math:`m_v` is the mass of water vapor, :math:`m_c` is the mass of liquid water, and :math:`m_p` is the mass of precipitate.
+From the definitions of the mass mixing ratio (ratio of mass of a component to mass of dry air), we have for any component
+
+.. math::
+    q_i \equiv = \frac{m_i}{m_a}.
+
+Using this we can write
+
+.. math::
+    \rho = m_a\frac{(1 + q_v + q_c + q_p)}{V}
+          = \rho_d(1 + q_v + q_c + q_p),
+
+where :math:`\rho_d \equiv \cfrac{m_a}{V}` is the density of dry air.
+
+
 Buoyancy
 =========
 
diff --git a/Source/TimeIntegration/ERF_make_buoyancy.cpp b/Source/TimeIntegration/ERF_make_buoyancy.cpp
index 7ad36127c..6e3b7f39e 100644
--- a/Source/TimeIntegration/ERF_make_buoyancy.cpp
+++ b/Source/TimeIntegration/ERF_make_buoyancy.cpp
@@ -173,9 +173,9 @@ void make_buoyancy (Vector<MultiFab>& S_data,
             amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
             {
                 Real rhop_hi = cell_data(i,j,k  ,Rho_comp) * (1.0 + qv_data(i,j,k  ) + qc_data(i,j,k  )
-                                                                  + qi_data(i,j,k  )) - r0_arr(i,j,k  );
+                                                                  + qi_data(i,j,k  )) + cell_data(i,j,k,RhoQp_comp) - r0_arr(i,j,k  );
                 Real rhop_lo = cell_data(i,j,k-1,Rho_comp) * (1.0 + qv_data(i,j,k-1) + qc_data(i,j,k-1)
-                                                                  + qi_data(i,j,k-1)) - r0_arr(i,j,k-1);
+                                                                  + qi_data(i,j,k-1)) +  cell_data(i,j,k-1,RhoQp_comp) - r0_arr(i,j,k-1);
                 buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo );
             });
         } // mfi