  An improved description for nonlinear plasma wakefields with phase velocities
near the speed of light is presented and compared against fully kinetic
particle-in-cell simulations. These wakefields are excited by intense particle
beams or lasers pushing plasma electrons radially outward, creating an ion
bubble surrounded by a sheath of electrons characterized by the source term $S
\equiv -\frac{1}{en_p}(\rho-J_z/c)$ where $\rho$ and $J_z$ are the charge and
axial current densities. Previously, the sheath source term was described
phenomenologically with a positive-definite function, resulting in a positive
definite wake potential. In reality, the wake potential is negative at the rear
of the ion column which is important for self-injection and accurate beam
loading models. To account for this, we introduce a multi-sheath model in which
the source term, $S$, of the plasma wake can be negative in regions outside the
ion bubble. Using this model, we obtain a new expression for the wake potential
and a modified differential equation for the bubble radius. Numerical results
obtained from these equations are validated against particle-in-cell
simulations for unloaded and loaded wakes. The new model provides accurate
predictions of the shape and duration of trailing bunch current profiles that
flatten plasma wakefields. It is also used to design a trailing bunch for a
desired longitudinally varying loaded wakefield. We present beam loading
results for laser wakefields and discuss how the model can be improved for
laser drivers in future work. Finally, we discuss differences between the
predictions of the multi- and single-sheath models for beam loading.
