The $(d+1)$-dimensional KPZ equation is the canonical model for the growth of rough $d$-dimensional random surfaces. A deep mathematical understanding of the KPZ equation for $d=1$ has been achieved in recent years, and the case $d\ge 3$ has also seen some progress. The most physically relevant case of $d=2$, however, is not very well-understood mathematically, largely due to the renormalization that is required: in the language of renormalization group analysis, the $d=2$ case is neither ultraviolet superrenormalizable like the $d=1$ case nor infrared superrenormalizable like the $d\ge 3$ case. Moreover, unlike in $d=1$, the Cole-Hopf transform is not directly usable in $d=2$ because solutions to the multiplicative stochastic heat equation are distributions rather than functions. In this article we show the existence of subsequential scaling limits as $\varepsilon \to 0$ of Cole-Hopf solutions of the $(2+1)$-dimensional KPZ equation with white noise mollified to spatial scale $\varepsilon$ and nonlinearity multiplied by the vanishing factor $|\log\varepsilon|^{-1/2}$. We also show that the scaling limits obtained in this way do not coincide with solutions to the linearized equation, meaning that the nonlinearity has a non-vanishing effect. We thus propose our scaling limit as a notion of KPZ evolution in $2+1$ dimensions.