Needs["NDSolve`FEM`"]

(*constants*)
ballRadiusSquared = 0.5;
rmax = 10.;
zmin = 0.;
zinit = 4.;
pa = 1.;
g = 9.8;
rho = 1.;
tmax = 0.2;

(*time dependent variable*)
zBallCentre[t_] := -0.5*g*t^2 + zinit + Sqrt[ballRadiusSquared];

(*define region coordinates*)
n = 1000;
(*circle coordinates*)
circTbl = Table[{Sqrt[ballRadiusSquared]*Cos[t], Sqrt[ballRadiusSquared]*Sin[t]}, {t, -Pi/2, Pi/2 - 2 Pi/(2 n), 2 Pi/(2 n)}];
Do[circTbl = ReplacePart[circTbl, {{i, 2}} -> circTbl[[i, 2]] + zBallCentre[0.7]], {i, Length[circTbl]}];
Do[If[circTbl[[i, 2]] > zinit, circTbl = Delete[circTbl, i], circTbl = circTbl], {i, Length[circTbl], 1, -1}];

(*surface coorinates*)
 surfaceTbl = Table[{x, zinit}, {x, 0., rmax, 1/n}];
If[Length[circTbl] == n, surfaceTbl = surfaceTbl,
    Do[If[circTbl[[Length[circTbl], 1]] > surfaceTbl[[1, 1]],surfaceTbl = Delete[surfaceTbl, 1],surfaceTbl = surfaceTbl],
    {i, 1, n}]];

(*put everything together*)
coordTbl = Join[circTbl, surfaceTbl];
coordTbl = Prepend[coordTbl, {0., zmin}];
coordTbl = Append[coordTbl, {rmax, zmin}];
lineTbl = Partition[Range[Length[coordTbl]], 2, 1];
lineTbl = Prepend[lineTbl, {Length[coordTbl], 1}];

(*mesh the region where de PDE will be solved*)
bmesh = ToBoundaryMesh["Coordinates" -> coordTbl, "BoundaryElements" -> {LineElement[lineTbl]}];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 0.1];
mesh["Wireframe"]

DirCond = DirichletCondition[[Phi][r, z] == 0., r == rmax];
NeumCond = NeumannValue[0., r == 0] + NeumannValue[0., r == rmax] + NeumannValue[0., z == zmin];

solution = NDSolveValue[{D[r*D[[Phi][r, z], r], r]/r + D[[Phi][r, z], {z, 2}] == NeumCond, DirCond}, [Phi], {r, z} [Element] mesh, "ExtrapolationHandler" -> {Automatic, "WarningMessage" -> False}]