We study the soft excess variability of the narrow line Seyfert 1 galaxy IRAS 13224-3809. We considered all five archival XMM-Newton observations, and we applied the flux-flux plot (FFP) method. We found that the flux-flux plots were highly affected by the choice of the light curves time bin size, most probably because of the fast and large amplitude variations, and the intrinsic non-linear flux--flux relations in this source. Therefore, we recommend that the smallest bin-size should be used in such cases. Hence, We constructed FFPs in 11 energy bands below 1.7 keV, and we considered the 1.7-3 keV band, as being representative of the primary emission. The FFPs are reasonably well fitted by a power-law plus a constant model. We detected significant positive constants in three out of five observations. The best-fit slopes are flatter than unity at energies below $sim 0.9$ keV, where the soft excess is strongest. This suggests the presence of intrinsic spectral variability. A power-law-like primary component, which is variable in flux and spectral slope (as $Gammapropto N_{rm PL}^{0.1}$) and a soft-excess component, which varies with the primary continuum (as $F_{rm excess}propto F_{rm primary}^{0.46}$), can broadly explain the FFPs. In fact, this can create positive `constants, even when a stable spectral component does not exist. Nevertheless, the possibility of a stable, soft--band constant component cannot be ruled out, but its contribution to the observed 0.2-1 keV band flux should be less than $sim 15$ %. The model constants in the FFPs were consistent with zero in one observation, and negative at energies below 1 keV in another. It is hard to explain these results in the context of any spectral variability scenario, but they may signify the presence of a variable, warm absorber in the source.