We present a detailed X-ray timing analysis of the highly variable NLS1 galaxy, IRAS 13224-3809. The source was recently monitored for 1.5 Ms with XMM-Newton which, combined with 500 ks archival data, makes this the best studied NLS1 galaxy in X-rays to date. We apply standard time-and Fourier-domain techniques in order to understand the underlying variability process. The source flux is not distributed lognormally, as expected for all types of accreting sources. The first non-linear rms-flux relation for any accreting source in any waveband is found, with rms ∝ flux 2/3 . The light curves exhibit significant strong non-stationarity, in addition to that caused by the rms-flux relation, and are fractionally more variable at lower source flux. The power spectrum is estimated down to ∼ 10 −7 Hz and consists of multiple peaked components: a lowfrequency break at ∼ 10 −5 Hz, with slope α < 1 down to low frequencies; an additional component breaking at ∼ 10 −3 Hz. Using the high-frequency break we estimate the black hole mass M BH = [0.5 − 2] × 10 6 M , and mass accretion rate in Eddington units, m Edd ∼ > 1. The broadband PSD and accretion rate make IRAS 13224-3809 a likely analogue of Very-high/Intermediate state black hole X-ray binaries. The non-stationarity is manifest in the PSD with the normalisation of the peaked components increasing with decreasing source flux, as well as the low-frequency peak moving to higher frequencies. We also detect a narrow coherent feature in the soft band PSD at 7×10 −4 Hz, modelled with a Lorentzian the feature has Q ∼ 8 and an rms ∼ 3 %. We discuss the implication of these results for accretion of matter onto black holes.