  In 1989 Mason and Newman proved that there is a 1-1-correspondence between
self-dual metrics satisfying Einstein vacuum equation (in complex case or in
neutral signature) and pairs of commuting parameter depending vector fields
$X_1(\lambda),X_2(\lambda)$ which are divergence free with respect to some
volume form. Earlier (in 1975) Pleba\'nski showed instances of such vector
fields depending of one function of four variables satisfying the so-called I
or II Pleba\'nski heavenly PDEs. Other PDEs leading to Mason--Newman vector
fields are also known in the literature: Husain--Park (1992--94), Schief
(1996). In this paper we discuss these matters in the context of the web
theory, i.e. theory of collections of foliations on a manifold, understood from
the point of view of Nijenhuis operators. In particular we show how to apply
this theory for constructing new ``heavenly'' PDEs based on different normal
forms of Nijenhuis operators in 4D, which are integrable similarly to their
predecessors. Relation with the Hirota dispersionless systems of PDEs and the
corresponding Veronese webs, which was recently observed by
Konopelchenko--Schief--Szereszewski, is established in all the cases. We also
discuss some higher dimensional generalizations of the ``heavelny'' PDEs and
the existence of related vacuum Einstein metrics in 4D-case.
