I'm trying to solve very simple differential equation with perdiodic boundary conditions and non-constant initial condition, however the results I'm getting don't seem to be correct. Here's my code:

size = 10 (* size of the domain *)
tMax = 20000000 (* time of the experiment *) 
pde = D[h[x, y, t], t] == -0.7 
BConds = {h[0, y, t] == h[size, y, t], h[x, 0, t] == h[x, size, t]} (* periodic boundary conditions *) 
topography[x1_, x2_] = Piecewise[{{400./size*x1 - 100., size/4. < x1 < size/2.}, {300. - 400./size*x1, size/2. < x1 <
    3.*size/4.}, {100., x1 == size/2.}, {0, x1 <= size/4.}, {0, x2 >= 3.*size/4.}}] 
ICond = h[x, y, 0] == topography[x, y] (* initial condition as a function of x *) 
sol = NDSolve[Join[{pde}, BConds, {ICond}], h, {x, 0, size}, {y, 0, size}, {t, 0, tMax}] 
fun[x_, y_, t_] := h[x, y, t] /. sol[[1]]

I wound expect, because of how the initial condition is specified, the solution to be symmetrical with respect to size/2. That's true for e.g. x=size/4 and x=3size/4, meaning .fun[size/4,0,0]=fun[3size/4,0,0].

However that's not true for x=size/3 and x=2*size/3, for which the solution in my opinion also should have the same values, but it doesn't:

fun[size/3, 0, 0]
30.8642
fun[2*size/3, 0, 0]
31.6872

Why the values aren't the same?

I tried to use

Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> "FiniteElement"}}

and even though it gives results much closer to what I expect, it doesn't allow me to use periodic boundary conditions.

The other problem is when I try to plot the solution:

Manipulate[ContourPlot[fun[x, y, t], {x, 0, 10}, {y, 0, 10}, ColorFunction -> (ColorData["Rainbow"][1 #] &), PlotLegends -> Automatic], {t, 0, tMax}]

The plot, for t=0, displays regions of fun[x,y,t]<0 even though the minimum value for the initial condition is 0. I would guess it's because of how the solution is interpolated, but could someone explain it to me in more details?