The two-dimensional self-consistent impedance method is used to calculate the electromagnetic surface impedance above subsurface structures at very low frequencies. The method was derived from Faraday's and Ampere's Laws and results in a linear matrix equation where the right hand side of the equation corresponds to the source field introduced into the model as a fixed magnetic value. An air layer above the earth's surface is included to allow the scattered magnetic field to be calculated at the surface. The source field is applied above the earth's surface as a Dirichlet boundary condition, and a Neumann boundary condition is applied to all other boundaries in the solution space. The left hand side of the linear equation corresponds to the impedance matrix determined by discretising the solution space into two-dimensional rectangular pixels or cells bounded by lumped impedance elements, with values determined by the electromagnetic properties of the local media and the size of the pixel in the model. The resulting sparse matrix offers the flexibility of cells of any shape or size. Due to the large matrix dimensions, an iterative solver with a preconditioning technique was used to improve the speed, size and convergence of the solution. The efficient forward modelling has been applied to the analysis of a coal seam with various structural anomalies and line of oxidation along a line defined by 500 m with 0.5 m resolution. This improved technique allows in-field inverse modelling of surface impedance data. This paper reports several likely coal-seam scenarios relevant to surface mining operations.