After the introduction of the ionization-injection scheme in laser wake field acceleration and of related high-quality electron beam generation methods, such as two-color and resonant multi-pulse ionization injection (ReMPI), the theory of thermal emittance has been used to predict the beam normalized emittance obtainable with those schemes. We recast and extend such a theory, including both higher order terms in the polynomial laser field expansion and non-polynomial corrections due to the onset of saturation effects on a single cycle. Also, a very accurate model for predicting the cycle-averaged distribution of the extracted electrons, including saturation and multi-process events, is proposed and tested. We show that our theory is very accurate for the selected processes of ${\mathrm{Kr}}^{8^{+}\to {10}^{+}}$ and ${\mathrm{Ar}}^{8^{+}\to {10}^{+}}$ , resulting in a maximum error below 1%, even in a deep-saturation regime. The accurate prediction of the beam phase-space can be implemented, for example, in laser-envelope or hybrid particle-in-cell (PIC)/fluid codes, to correctly mimic the cycle-averaged momentum distribution without the need for resolving the intra-cycle dynamics. We introduce further spatial averaging, obtaining expressions for the whole-beam emittance fitting with simulations in a saturated regime, too. Finally, a PIC simulation for a laser wakefield acceleration injector in the ReMPI configuration is discussed.