  Using a combinatorial method, the partition functions for two-dimensional
nearest neighbour Ising models have been derived for a square lattice of 16
sites in the presence of the magnetic field. A novel hierarchical method of
enumeration of all the configurations for any arrangement of sites has been
proposed. This enumeration has been executed by a systematic analysis of the
appropriate diagrams without employing any algorithmic approach or
computational tools. The resulting algebraic eqn in terms of the magnetic field
and nearest neighbour interaction energies may then provide a methodology for
deducing the magnetization in the thermodynamic limit of infinite sites. A
semi-empirical eqn for magnetization is proposed for non-zero magnetic fields.
