Calculation for equation (4.1), which gives the local factor for k=3, l > 0, and any composite shift h. (For the codes below, h is \nu_p(h).)

c3[ell_,h_]:=Simplify[(1+(2*Sum[Sum[Binomial[j+j2+ell-1,ell-1]/p^(j),{j,0,Infinity}]/p^(j2),{j2,1,h}]*(1-1/p)^ell+Sum[Sum[Sum[Binomial[j+j1+j2+ell-1,ell-1]/p^(j),{j,0,Infinity}]/p^(j2),{j2,1,h-j1}]/p^j1,{j1,1,h-1}]*(1-1/p)^ell+Binomial[h+ell-1,ell-1]*(1-1/p)^(ell-1)*Sum[Sum[1/p^(j2),{j2,0,Infinity}]/p^(j1),{j1,h+1,Infinity}]+Binomial[h+ell-1,ell-1]*(1-1/p)^(ell-1)*Sum[Sum[1/p^(j2*(1+0)),{j2,h+1-j1,Infinity}]/p^(j1*(0+1)),{j1,0,h}]))/( 1 + (1-1/p)^(ell-1)*Sum[Binomial[2,j]/(p-1)^j,{j,1,2}] )];
c3[3,h]

c3[4,h]

Compare with the predicted local factor from, e.g., equation (1.7) of Ng and Thom:

f[k_,ell_,h_]:=Simplify[( 1 + Sum[(Binomial[i+k-1,k-1]*Binomial[i+ell-1,ell-1]-Binomial[i+k-2,k-1]*Binomial[i+ell-2,ell-1])/p^i,{i,1,h}]  + Sum[(Binomial[h+k-1,k-1]*Binomial[i+ell-2,ell-2]+ Binomial[h+ell-1,ell-1]*Binomial[i+k-2,k-2])/p^i,{i,h+1,Infinity}] )/((1-1/p)^(-k+1)+(1-1/p)^(-ell+1)-1)];

f[3,3,h]

f[3,4,h]

The difference c3[ell,h] -  f[3,ell,h] for ell<=5, showing that the two constants match:

Table[c3[ell,h]-f[3,ell,h],{ell,5}]