We simulate the multiphase interstellar medium (ISM) randomly heated and stirred by supernovae (SNe), with gravity, differential rotation and other parameters of the solar neighbourhood. Here we describe in detail both numerical and physical aspects of the model, including injection of thermal and kinetic energy by SN explosions, radiative cooling, photoelectric heating and various transport processes. With a three-dimensional domain extending 1 × 1 kpc2 horizontally and 2 kpc vertically (symmetric about the galactic mid-plane), the model routinely spans gas number densities 10−5-102 cm−3, temperatures 10-108 K and local velocities up to 103 km s−1 (with Mach number up to 25). The working numerical resolution of 4 pc has been selected via simulations of a single expanding SN remnant, where we closely reproduce, at this resolution, analytical solutions for the adiabatic and snowplough regimes. The feedback of the halo on the disc cannot be captured in our model where the domain only extends to the height of 1 kpc above the mid-plane. We argue that to reliably model the disc-halo connections would require extending the domain horizontally as well as vertically due to the increasing horizontal scale of the gas flows with height. The thermal structure of the modelled ISM is classified by inspection of the joint probability density of the gas number density and temperature. We confirm that most of the complexity can be captured in terms of just three phases, separated by temperature borderlines at about 103 and 5 × 105 K. The probability distribution of gas density within each phase is approximately lognormal. We clarify the connection between the fractional volume of a phase and its various proxies, and derive an exact relation between the fractional volume and the filling factors defined in terms of the volume and probabilistic averages. These results are discussed in both observational and computational contexts. The correlation scale of the random flows is calculated from the velocity autocorrelation function; it is of the order of 100 pc and tends to grow with distance from the mid-plane. We use two distinct parametrizations of radiative cooling to show that the multiphase structure of the gas is robust, as it does not depend significantly on this choice.