The discovery of a highly-charged sheath region at the boundary between a plasma and a surface is one of the earliest and most important discoveries in plasma science. However sheath physics has almost always been omitted from studies of the dynamics of plasma-facing liquid surfaces which are rapidly assuming a pivotal role in numerous industrial and fusion applications. This paper presents full simulations of the plasma-sheath-liquid interface and finds good agreement with theoretical stability limits and experimental observations of cone formation and pulsed droplet ejection. Consideration of sheath physics is strongly encouraged in all future studies of plasma-liquid interactions.Plasma-liquid interactions occur in a diverse and ever-increasing range of technological applications such as advanced materials processing, spacecraft propulsion, nanoparticle synthesis, chemical catalysis, plasma medicine, food treatment and water purification [1-3]. Furthermore there is increasing recognition of their role in magnetic fusion devices where the intense heat loads of the fusion plasma can cause melting and permanent damage to the metal walls [4,5]. One intriguing development is to coat the plasma-facing walls with liquid lithium which enables self-replenishment of any eroded material [6][7][8]. The production of liquid droplets is observed from both accidentally-melted and intentionally-molten metal surfaces [9,10].Maintaining the stability of the plasma-liquid interface remains a crucial challenge in these new technologies and, accordingly, has received widespread attention: Kelvin-Helmholtz [11], Rayleigh-Taylor [12] and Rayleigh-Plateau [13] instabilities, as well as droplet ejection from liquid splashing [14] or the bursting of surface bubbles [15], have all been identified as potential causes of surface disruption and droplet ejection. However all of these studies have neglected the fundamental role of the plasma sheath at the plasma-liquid interface. Langmuir's discovery of this sheath region, and the large electric fields contained within it, is one of the earliest discoveries in plasma science [16]; it is, therefore, surprising that sheath effects have only recently been theoretically investigated with a linear perturbation analysis [17]. This paper expands on this theoretical work by developing complementary numerical simulations. Comparisons are drawn between theory, simulation and relevant experimental observations.The linear perturbation theory in [17] took the simple Bohm model of a plasma sheath [18,19] which comprises a collisionless cold-ion fluid, with density n, velocity u and individual ion mass m i , and Boltzmanndistributed electrons, with temperature T e and sheath-edge density n 0 , which interact with each other and a conducting wall via the electrostatic field E as described by the potential f. The corresponding sheath equations are * Substantial portions of this paper are adapted from section 4.2 of the first author's recent PhD thesis.