Boolean networks are simplified models of gene regulatory networks. We derive an approximation of the size distribution of perturbation avalanches in Boolean networks based on known results in the theory of branching processes. We show numerically that the approximation works well for different kinds of Boolean networks. It has been suggested that gene regulatory networks may be dynamically critical. To study this, as an application of the presented theory we present a novel method for estimating an order parameter from microarray data. According to the available data and our method, we find that gene regulatory networks appear to be stable and reside near the phase transition between order and chaos.