Energy storage is \$E = \frac{1}{2} L I^2\$, when the system is linear (L is constant). Otherwise, we must make it a function of I and integrate both.
Start with the fundamentals.
First, note that inductance is the incremental value, meaning the inductance measured for small changes in current and flux. That is,
$$ L \equiv \frac{d\phi}{di} $$
Energy is defined as,
$$ E = \int v(t) \, i(t) \, dt $$
for given voltage and current waveforms. We can test with a simple ramp,
$$ i(t) = \frac{I_0 t}{t_0} $$
and integrate time from 0 to t0. We then need the voltage. Note \$d\phi = v(t) dt\$, so we have
$$ v(t) = \frac{L(i(t)) di}{dt} $$
Substitute in i(t), and its derivative (which is trivial), to get L and thus v as functions of time.
We need a function for L, to symbolically evaluate this further. Suppose we choose a simple polynomial sigmoid function (Lorentzian),
$$ L(i) = \frac{a}{i^2 + b^2} + L_0 $$
where \$\frac{a}{b^2} + L_0\$ is the zero-bias inductance, \$L_0\$ is the fully-saturated inductance, and \$b\$ determines the saturation point (the half-way point between nominal and minimal permeability is at \$i = b\$). b has units of A, and a is H A-2.
Substituting,
$$ E = \int_0^{t_0} L\left(\tfrac{I_0 t}{t_0}\right) \frac{I_0}{t_0} \, \frac{I_0 t}{t_0} \, dt \\
E = \frac{{I_0}^2}{{t_0}^2} \int_0^{t_0} \left( \frac{a t}{(\frac{I_0 t}{t_0})^2 + b^2} + L_0 t\right) \, dt \\
E = \frac{{I_0}^2}{{t_0}^2} \int_0^{t_0} \left( \frac{a \, t \, {t_0}^2}{{I_0}^2 t^2 + b^2 {t_0}^2} + L_0 t\right) \, dt \\
E = \frac{{I_0}^2 a \, {t_0}^2}{{I_0}^2 {t_0}^2} \int_0^{t_0} \frac{t}{t^2 + \frac{b^2 {t_0}^2}{{I_0}^2}} \, dt + \frac{{I_0}^2 L_0}{{t_0}^2} \int_0^{t_0} t \, dt \\
E = \frac{a}{2} \ln \left( \frac{{I_0}^2}{b^2} + 1 \right) + \frac{L_0 {I_0}^2}{2} $$
What does this mean? Well, we can approximate \$\ln(x+1)\$ with the Maclaurin series \$x + \frac{x^2}{2} - \frac{x^3}{3} + \cdots\$; in particular, at small currents (i.e., truncating to the first term), it reduces to
$$ E = \frac{{I_0}^2}{2} \left( \frac{a}{b^2} + L_0 \right) $$
which tracks with our initial inductance value.
Asymptotically, \$\ln\$ is a slowly-increasing function, and \$x^2\$ quickly outpaces it, so the \$L_0\$ term dominates in the long run -- as we expect for a saturated inductor.
The saturation is a bit hard to see with this formula; this is more or less as expected, because a Lorentzian has a rather subtle drop-off (even powdered iron inductors have, I think, a steeper curve than this). It's further complicated by the energy curve itself (it's always upwards-trending, we can only look for inflection points in it). But we still observe something reasonable when we plot the function. Consider some basic values of the form:
Source: https://www.wolframalpha.com/input?i=50+ln%28x%5E2%2B1%29%2Bx%5E2
At first, energy rises relatively quickly with current, in the unsaturated region (lots of flux); eventually, the rate subsides as saturation sets in. At large values, it keeps on going, despite the low inductance, because current is just so high up there.
Finally, to address the question:
Energy storage will always be increasing; at least, as long as there isn't any mechanical movement to sap energy from the system. The quadratic term must dominate, in the end.
The real question is this:
How much current can a circuit reasonably use?
The more general take-away is, I think, in terms of ratios, of power matching. Consider this: how do we determine how many turns to put on the core? We have a core of some nominal energy storage capacity, independent of the winding we put around it; the question, then, is what impedance -- what ratio of voltage to current -- the circuit needs. The turns count is the transformer matching ratio for an inductor to the circuit.
And when the inductor saturates, that ratio changes. A suddenly-low inductance is no longer useful to a circuit designed around a given impedance; less energy can be stored up to a given current limit, or within the given time scale the circuit operates at (Fsw).
There are further extenuating factors that restrict the design range of an inductor. An electronic filter needs to land more-or-less on target frequency; a control system needs to remain stable and responsive. There are some circumstances where a varying inductor can be useful (e.g. "swinging choke" in an otherwise-passive power supply, no control issues to worry about), but mostly we merely tolerate the variation, and desire a modestly bounded range, perhaps as wide as 50-150% for SMPS control and filtering purposes, or much tighter for various filter applications (a precision radio filter might need much better than ±1%).
