An exactly solvable model of 2D Majorana zero modes

Well it depends on what you mean by "exactly solvable". Many people would consider any quadratic Hamiltonian(i.e. non-interacting) exactly solvable, so the BCS Hamiltonian of a p+ip superconductor already falls into this category. One can demand more -- like at a special point of Kitaev chain ($t=\Delta, \mu=0$), the Hamiltonian consists of commuting terms. This is in my opinion the strongest form of exact solubility. However we believe that this is not possible for chiral p+ip superconductors. There are other exactly solvable models (in the sense of commuting terms), in which certain types of vortices/fluxes bind Majorana zero modes.