Domain structures of small pseudo-single domain (PSD) magnetite near the Verwey transition (Tυ) at ≈120 K were modeled using an unconstrained three-dimensional micromagnetic algorithm. The single-domain (SD) threshold (d0) for the monoclinic phase below Tυ was calculated to be ≈0.14 μm at 110 K. However, it is postulated that as a result of the very high energy barriers in the monoclinic phase, grains near d0 in size and in vortex states are unlikely to denucleate domain walls to become SD. Low-temperature cycling of saturation isothermal remanence (SIRM), thermoremanence (TRM), and partial TRM (pTRM) through Tυ was simulated. Domain structures were found to align along the monoclinic “easy” magnetocrystalline anisotropy axis, i.e., the c axis, on simulated cooling through Tυ. This process was found to “destroy” SIRM structures giving rise to demagnetization; however, for TRM and pTRM structures only “closure” domains were removed increasing magnetostatic leakage giving rise to a reversible anomaly in rough agreement with experimental studies. SIRM displayed a smaller anomaly at Tυ, in agreement with experimental studies.