Willis fluids are characterized by constitutive relations that couple the pressure and momentum density to both the particle velocity and the volume strain. This effective dynamic response coupling may arise due to microstructural asymmetry, long range order, or time-varying material properties and has been shown to be analogous to electromagnetic bianisotropy [Phys. Rev. B 96, 104303 (2017)]. Recent work has shown that, when wavevectors exceed a critical magnitude, the currently employed model for passive Willis fluids possesses a linear instability stemming from the truncation of higher-order coupling terms in the derivation of the constitutive relationships [J. Acoust. Soc. Am. 144(3), 1832 (2018)]. While this instability can be avoided when using frequency domain methods, it poses significant problems for time domain simulations, as short-wavelength numerical noise may grow exponentially. In the present work, we report on methods of stabilizing the Willis fluid model. Specifically, we consider mathematical manipulations of the wave equation derived from the lowest-order constitutive equations, similar to methods that have been explored in continuum approximations of nonlinear spring-mass chains, as well as higher-order coupling terms. Physical interpretations and implications for time domain finite difference and finite element simulations are discussed. [Work supported by the postdoctoral program at ARL:UT.]Willis fluids are characterized by constitutive relations that couple the pressure and momentum density to both the particle velocity and the volume strain. This effective dynamic response coupling may arise due to microstructural asymmetry, long range order, or time-varying material properties and has been shown to be analogous to electromagnetic bianisotropy [Phys. Rev. B 96, 104303 (2017)]. Recent work has shown that, when wavevectors exceed a critical magnitude, the currently employed model for passive Willis fluids possesses a linear instability stemming from the truncation of higher-order coupling terms in the derivation of the constitutive relationships [J. Acoust. Soc. Am. 144(3), 1832 (2018)]. While this instability can be avoided when using frequency domain methods, it poses significant problems for time domain simulations, as short-wavelength numerical noise may grow exponentially. In the present work, we report on methods of stabilizing the Willis fluid model. Specifically, we consider mathema...