Color of line in ParametricPlot according to indepedent variable

I solve a set of pdes numerically using NDSolve:

ic1[x_, y_] := Which[x < 3 && y < 3, 6, True, 0] ic2[x_, y_] := Which[x < 3 && y < 3, 0.5, True, 1/c11] sol = NDSolve[{D[pp[t, x, y], t] ==  0.05*Laplacian[pp[t, x, y], {x, y}] +   pp[t, x,     y]*(1 - c11*pp[t, x, y] - z[t, x, y]/(1 + pp[t, x, y]^2)),  D[z[t, x, y], t] ==  0.05*Laplacian[z[t, x, y], {x, y}] +   z[t, x, y]*(eps*pp[t, x, y]/(1 + pp[t, x, y]^2) - m), (D[    pp[t, x, y], x] /. x -> 0) ==  0, (D[pp[t, x, y], y] /. y -> 0) ==  0, (D[z[t, x, y], x] /. x -> 0) ==  0, (D[z[t, x, y], y] /. y -> 0) ==  0, (D[pp[t, x, y], x] /. x -> 20) ==  0, (D[pp[t, x, y], y] /. y -> 20) ==  0, (D[z[t, x, y], x] /. x -> 20) ==  0, (D[z[t, x, y], y] /. y -> 20) == 0, z[0, x, y] == ic1[x, y],  pp[0, x, y] == ic2[x, y]}, {pp, z}, {t, 0,  100} {x, 0, 20}, {y, 0, 20}]; 

First of all, evaluation over a long time interval does not work. Can I do something to export only certain time steps to save memory or something?

Furthermore, I want to plot the results of a certain point in space (x,y) in a parametric plot. I can do this with the following code:

ParametricPlot[ Evaluate[{pp[t, 4, 4], z[t, 4, 4]} /. sol], {t, 0, 100},  PlotRange -> {{0, 2}, {0, 4}}, AxesLabel -> {p, z},  PlotStyle -> {Thickness[0.002]},  ColorFunction -> {Function[{x, y, t}, Hue[t]]}] 

However, I find it hard to link the Hue-Colors to the temporal evolution. I would prefer a color gradient with changing opacity or some temperature (CoolWarm) scale or something. How can I do this? I thought this should be quite straight forward but I didn’t find anything about it the last days.

Thank you!