Physics Department, University of California, Davis, CA, 95616(Dated: July 15, 2011)We compute the bipartite entanglement properties of the spin-half square-lattice Heisenberg modelby a variety of numerical techniques that include valence bond quantum Monte Carlo (QMC),stochastic series expansion QMC, high temperature series expansions and zero temperature couplingconstant expansions around the Ising limit. We find that the area law is always satisfied, but inaddition to the entanglement entropy per unit boundary length, there are other terms that dependlogarithmically on the subregion size, arising from broken symmetry in the bulk and from theexistence of corners at the boundary. We find that the numerical results are anomalous in severalways. First, the bulk term arising from broken symmetry deviates from an exact calculation that canbe done for a mean-field N´eel state. Second, the corner logs do not agree with the known results fornon-interacting Boson modes. And, third, even the finite temperature mutual information shows ananomalous behavior as T goes to zero, suggesting that T → 0 and L → ∞ limits do not commute.These calculations show that entanglement entropy demonstrates a very rich behavior in d > 1,which deserves further attention.I. INTRODUCTION