Discrepancy between the results of NIntegrate with different methods and options

I am trying to perform a numerical integration on a function defined through a sum of exponential terms. The summation is given by:

sum[z_, z0_, t_, nmax_] := 
 1/Sqrt[4 t]*Sum[ Exp[-(z - z0 - 2 n)^2/(4 t)] + Exp[-(z + z0 - 2 n)^2/(4 t) ], {n, -nmax, nmax}];

and we define

f[z_, zp_, yp_, z0_, t_] = 
  ( Exp[-(yp)^2/(8t)]/ (8t) )*D[sum[z, z0, t, 10]*sum[zp, z0, t, 10], t, z, zp];

where I have chosen nmax=10 (update I tried a smaller value for nmax and still got similar behaviour. In principle, nmax should be infinity for my calculation, but in that case the sum evaluates to EllipticTheta[3,,] and integration with that has its own issues as in this post)

I wish to perform numerical integration on f. I define

int1[a_?NumericQ] := 
 NIntegrate[ f[zp, zpp, ypp, z0, t]*( ypp (zp - zpp) )/( a (zp - zpp)^2 + ypp^2 )^(3/2), 
{z0, 0., 1.}, {t, 0., 10.}, {zp, 0., 1.}, {zpp, 0., 1.}, {ypp, 0., Infinity}]

Based on which I get the following table (no errors generated – update: after restarting the kernel, it did generate NIntegrate::slowcon errors)

tab1 = Table[{a, int1[a]}, {a, 1., 5., .5}]

(* {{1., 0.0577262}, {1.5, -0.0245604}, {2., -0.000784569},
{2.5,0.000552247}, {3., -246.874}, {3.5, -169.718}, 
{4., -118.106}, {4.5,0.0452358}, {5., -0.503675}} *)

To check, I tried AccuracyGoal-> 30 and PrecisionGoal -> 30 and got the same results. However, as soon as I include WorkingPrecision, the integrated gives 0 all the time. For example:

intwrk[a_] := 
 NIntegrate[ f[zp, zpp, ypp, z0, t]*(ypp (zp - zpp))/(a (zp - zpp)^2 + ypp^2)^(3/2), 
{z0, 0, 1}, {t, 0, 10}, {zp, 0, 1}, {zpp, 0, 1}, {ypp, 0, Infinity}, 
AccuracyGoal -> 30, PrecisionGoal -> 30, WorkingPrecision -> 100]

Gives all zeros – which I can’t seem to understand. Am I doing something wrong here?

(update: there is now a NIntegrate::precw error stating the precision of the argument is less than the working precision)

Then I also performed the calculation using LocalAdaptive method:

int2[a_] := 
 NIntegrate[ f[zp, zpp, ypp, z0, t]*(ypp (zp - zpp))/(a (zp - zpp)^2 + ypp^2)^(3/2),
{z0, 0., 1.}, {t,0., 10.}, {zp, 0, 1}, {zpp, 0., 1.}, {ypp, 0., Infinity }, Method -> "LocalAdaptive"]

which gives

tab2 = Table[{a, int2[a]}, {a, 1., 5., .5}]

(* {{1., -2.87734*10^-19}, {1.5, 8.24258*10^-20}, {2., -8.72407*10^-20},
{2.5, -2.78288*10^-20},{3., -3.98938*10^-20}, {3.5, -3.81874*10^-21},
{4., -2.05326*10^-20},{4.5, 8.15403*10^-21}, {5., -2.85668*10^-20}} *)

which is a very different result. This one also runs much faster than the int1. I guess these numbers here are not reliable, but how can I check them?

PS: based on the answer to my previous post, I tried Method -> "GaussKronrodRule" and Method -> {"MultidimensionalRule", "Generators" -> 9}, but they run forever and I couldn’t get any outcome.