We consider elliptic divergence form operators $L$ in open subsets of Euclidean space. We let the Dirichlet energy induce a heat flow simultaneously in the set and on the boundary. This corresponds to dynamic boundary conditions in the parabolic problem. We will explain how to prove $L^p$ square function estimates for the corresponding heat semigroup $\exp(tL)$. The main obstacle is that different parts of the underlying geometric domain (namely, interior and boundary) have different scaling dimensions.