memory - How to delete InterpolationFunction created by NDSolve/ParametricNDSolve?

I have faced a uncommon problem. The situation is that I need use ParametricNDSolve/ParametricNDSolveValue to solve ODEs system in the non-linear model fitting process. So every iteration generating a new parameter set and create a new InterpolatingFunction object. I guess the older InterpolatingFunction objects are not abandoned and it may be the reason for memory increase rapidly during the process. Now the problem is whether I can clear the objects at every iteration and how to do it.

For a similar problem, I copy a 4D PDE example from MathGroup Archive:

i = 0;

sol = NDSolve[{D[u[t, x, y, z], t] == 
     D[u[t, x, y, z], x, x] + D[u[t, x, y, z], y, y] + 
     D[u[t, x, y, z], z, z], u[0, x, y, z] == 0, 
     u[t, 0, y, z] == Sin[t], u[t, 40, y, z] == 0, 
     u[t, x, 0, z] == Sin[t], u[t, x, 40, z] == 0, 
     u[t, x, y, 0] == Sin[t], u[t, x, y, 40] == 0}, {u}, {t, 0, 
     100}, {x, 0, 40}, {y, 0, 40}, {z, 0, 40}, MaxSteps -> Infinity, 
    MaxStepSize -> 1, 
    EvaluationMonitor :> If[t > i, Print[{t, MemoryInUse[]/1024^2 // N}];
      i += 10;]]

MemoryInUse[]/1024^2 // N

This resulting InterpolatingFunction object will use about 1.5G memory. Clear[sol] is invalid to release its memory. If the PDE is used in model fitting, how to get the final result without memory running out?