  Accurate simulation of nonlinear acoustic waves is essential for the
continued development of a wide range of (high-intensity) focused ultrasound
applications. This article explores mixed finite element formulations of
classical strongly damped quasilinear models of ultrasonic wave propagation;
the Kuznetsov and Westervelt equations. Such formulations allow simultaneous
retrieval of the acoustic particle velocity and either the pressure or acoustic
velocity potential, thus characterizing the entire ultrasonic field at once.
Using non-standard energy analysis and a fixed-point technique, we establish
sufficient conditions for the well-posedness, stability, and optimal a priori
errors in the energy norm for the semi-discrete equations. For the Westervelt
equation, we also determine the conditions under which the error bounds can be
made uniform with respect to the involved strong dissipation parameter. A
byproduct of this analysis is the convergence rate for the inviscid (undamped)
Westervelt equation in mixed form. Additionally, we discuss convergence in the
$L^q(\Omega)$ norm for the involved scalar quantities, where $q$ depends on the
spatial dimension. Finally, computer experiments for the Raviart--Thomas (RT)
and Brezzi--Douglas--Marini (BDM) elements are performed to confirm the
theoretical findings.
