diff --git a/docs/source_docs/user_guide/inputs/solids_model.rst b/docs/source_docs/user_guide/inputs/solids_model.rst
index 7b57e3b5a9823a4217e649e907a787935881e7a5..7b2fc2991fd757cb2a52873e52fe851a3803389d 100644
--- a/docs/source_docs/user_guide/inputs/solids_model.rst
+++ b/docs/source_docs/user_guide/inputs/solids_model.rst
@@ -87,6 +87,10 @@ Enabling the DEM solver and specifying model options.
 | dem.damping_tang_fac.pw | Factor relating the tangential damping coefficient to the normal        | Real     |  0.5      |
 |                         | damping coefficient :: particle to wall collisions                      |          |           |
 +-------------------------+-------------------------------------------------------------------------+----------+-----------+
+| dem.implicit_drag       | Apply fluid-particle drag force by                                      | int      |  0        |
+|                         | :ref:`implicit velocity update<InputsDEMImplicitDrag>`.                 |          |           |
++-------------------------+-------------------------------------------------------------------------+----------+-----------+
+
 
 The following inputs use the DEM type names specified using the `dem.solve` input to define restitution coefficients and
 are proceeded with `dem.restitution_coeff`. These must be defined for all solid-solid and solid-wall combinations.
@@ -126,6 +130,40 @@ Below is an example for specifying the inputs for two DEM solids.
    dem.restitution_coeff.sand.wall =  0.85
    dem.restitution_coeff.char.wall =  0.89
 
+.. _InputsDEMImplicitDrag:
+
+Implicit drag
+~~~~~~~~~~~~~
+
+By default, the fluid-particle drag force is applied explicitly to particles,
+
+.. math::
+
+   u_p^{n^{\prime}+1} = u_p^{n^\prime} + dt_{\mathrm{DEM}}\left(\beta^n (u_f^{n+1} - u_p^{n}) + \mathcal{F}_p^{n^\prime}\right)/m_p
+
+where
+:math:`u_p` is the `p`-th particle velocity at the indicated time level,
+:math:`\beta^n` is the drag coefficient,
+:math:`u_f^{n+1}` is the fluid velocity,
+:math:`\mathcal{F}_p` are all non-drag forces (e.g., gravitational, buoyancy, collision etc.), and
+:math:`m_p` is particle mass.
+Terms exist at different time levels due in part to how the fluid and particle models are coupled
+and partly because particles typically sub-step in time. Specifically, a time step advances the fluid from
+:math:`t^n` to :math:`t^{n+1}`, then particles take multiple smaller time steps to traverse the same
+total time. In the above equation, prime markers differentiate quantities updated at each DEM sub-step.
+
+The updated particle velocity, :math:`u_p^{n^{\prime}+1}`, is substituted into the drag expression when
+`implicit drag` is enabled.
+
+.. math::
+
+   u_p^{n^{\prime}+1} = \frac{u_p^{n^\prime} + dt_{\mathrm{DEM}}\left( \beta^n u_f^{n+1} + \mathcal{F}_p^{n^\prime} \right)/m_p}{1 + dt_{\mathrm{DEM}} \beta^n / m_p}
+
+By using the updated velocity, the drag force computed for the fluid and particles is inconsistent,
+and interphase momentum is not conserved. However, implicit drag may be needed for numerical stability
+when the particle density is much less than the fluid (such as bubbles).
+
+
 
 PIC model settings
 ------------------