This study is a follow-up to the numerical simulations of the formation of Lisan seismites via shear instability of stratified fluid (known as Kelvin Helmholtz Instability – KHI). The study and the simulations, performed by Wetzler, Marco, Agnon and Heifetz (2005, 2010,2020), suggest that KHI is a plausible mechanism of earthquake-induced shear to perform various patterns of the deformations, observed both in the Lisan outcrops, and in the Dead Sea deep drill cores. In those simulations, it was assumed that the sediments behave as a Newtonian fluid. In reality, however, sediments exhibit viscoplastic rheology with nonlinear relations between the stress and the deformation rate tensors. For small-applied shear stress the sediments behave as solids and only after overcoming a certain threshold (yield stress), the sediments begin to flow. Furthermore, sediments tend to obey a shear thinning behavior, meaning their resistance to the applied stress (i.e. their effective viscosity) decreases with the increase of the deformation rate. In the current study, we first measured the viscoplastic rheological characteristics of carefully prepared samples of Dead Sea sediments (aragonite, detritus and their mixture), by using a rotary shear. We found substantial expected differences between the ductile detritus and the brittle aragonite. These data were implemented in the viscoplastic Herschel-Bulkley model, relating the stress and the deformation rate tensors. Next we nested the Herschel-Bulkley model in the state-of-the-art Basilisk software to solve the full Navier-Stokes fluid equation, using a high-resolution adaptive Cartesian mesh solver. The viscoplastic simulations of KHI dynamics provide more realistic and detailed structures of various forms of the observed seismites. The seismites' deformation intensity is related to the peak ground acceleration and the layer thickness via the bulk Reynolds and Richardson numbers, which are now (as opposed to Newtonian dynamics) time dependent.