Heat diffusion - Is this the correct approach to include Newmann boundary conditions?

Thank you for looking at this problem. Is this the correct approach to include neumann boundary conditions? With this solution temperature is not correct, and there´s no diffusion. The model seems unbalanced. However if the area is included in km2 instead of m2, temperature is ok and there is diffusion...

I´m trying to aplly the following 1D heat diffusion equation to a lake:

$$\frac{∂T}{∂t}= \frac{1}{A(z)}\ \frac{∂}{∂z}\ ((A(z)K(z)\frac{∂T}{∂z})$$

Where: k(z)= Eddy diffusion coefficient, 1.43x10-3 (m2/s), A(z)= lake area at depth z, m2

The surface boundary condition is the following:

$$ρcK(z,t)\frac{∂T}{∂z}= -Q$$

Where: ρ= water density, 1000 kg/m3 c= water specific heat, 4186 J/kgºC Q= radiation, w/m2 (or J.s-1.m-2). The lake bottom condition is the following:

$$\frac{∂T}{∂z}= 0$$

I´m considering the Crank-Nicolson scheme with centered differences in space and time and I have used Thomas Algorithm to solve the system of equations.

The equation discretization is included in the picture below. But to simplify the analysis I´m including here the equations for each coefficient (f = k(z); n is time and i is space):

$$α_{i}=\frac{1}{A_{i,n+1}+A_{i,n}}$$ $$s=\frac{∆t}{2(∆z)^{2}}$$

$$a_{i} = −α_{i}s(f_{i−1,n+1}+ f_{i,n+1})$$

$$b_{i} = 1+α_{i}s(f_{i−1,n+1}+ 2f_{i,n+1}+f_{i+1,n+1})$$

$$c_{i} = −α_{i}s(f_{i,n+1}+ f_{i+1,n+1})$$

$$d_{i,n} = α_{i}s(f_{i,n}+ f_{i-1,n})T_{i−1,n}+α_{i}s(f_{i+1,n}+ f_{i,n})T_{i+1,n}+[1−α_{i}s(f_{i−1,n}+ 2f_{i,n}+f_{i+1,n})]T_{i,n}$$

for the surface boundary:

$$b_{0}=1+α_{0}s(f_{−1,n+1}+ 2f_{0,n+1}+f_{1,n+1})$$

$$c_{0}=−α_{0}s(f_{-1,n+1}+ 2f_{0,n+1}+f_{1,n+1})$$

$$d_{0,n}=[1-α_{0}s(f_{-1,n}+ f_{0,n}+f_{1,n}]T_{0,n}+α_{0}s(f_{-1,n}+2 f_{0,n}+f_{1,n})T_{1,n}+α_{0}s(f_{−1,n}+ 2f_{0,n})\frac{2Q_{n}∆z}{ρcK_{0,n}}+α_{0}s(f_{-1,n+1}+ f_{0,n+1})\frac{2Q_{n+1}∆z}{ρcK_{0,n+1}}$$

for the bottom boundary:

$$a_{NZ} = −α_{NZ}s(f_{NZ−1,n+1}+ 2f_{NZ,n+1}+f_{NZ+1,n+1})$$

$$b_{NZ} = 1+α_{NZ}s(f_{NZ−1,n+1}+ 2f_{NZ,n+1}+f_{NZ+1,n+1})$$

$$d_{NZ,n} = α_{NZ}sT_{NZ−1,n}(f_{NZ-1,n}+ 2f_{NZ,n}+f_{NZ+1,n})+[1−α_{NZ}s(f_{NZ−1,n}+ 2f_{NZ,n}+f_{NZ+1,n})]T_{NZ,n}$$

I have also include a simplification of the Python code. I have removed all the loops and I have considered a constante value for eddy diffusivity K(z). It´s a daily simulation (N=10) for a lake with 6 layers. Each layer is 2m deep (z=2m).

import numpy as np
np.set_printoptions(threshold=np.inf)
import pylab as plt

#------------------------------------------------------------------------------------------------------------------------------------
# Tri Diagonal Matrix Algorithm (Thomas algorithm) solver (source:https://gist.github.com/cbellei/8ab3ab8551b8dfc8b081c518ccd9ada9)
#-------------------------------------------------------------------------------------------------------------------------------------

def TDMAsolver(a, b, c, d):

nf = len(d)
ac, bc, cc, dc = map(np.array, (a, b, c, d))
for it in xrange(1, nf):
mc = ac[it-1]/bc[it-1]
bc[it] = bc[it] - mc*cc[it-1]
dc[it] = dc[it] - mc*dc[it-1]

xc = bc
xc[-1] = dc[-1]/bc[-1]

for il in xrange(nf-2, -1, -1):
xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]

del bc, cc, dc

return xc

#-------------------------------------------------------------------------------------------------------------------------------------
#Crank Nicolson scheme
#-------------------------------------------------------------------------------------------------------------------------------------
#Simulation time period (10 days)
N=10

Nlayers=6

#Space, m
DZ=2.0

Q=300.0

#Eddy diffusion coefficient at 25 C (m2/s)
K=1.43*10**-3

#Water temperature array
TempLake=np.zeros((N+1,Nlayers))

# Initial condition, C
TempLake[0]=[25,25,25,25,25,25]

#Lake surface area, m2
A1=23992000.0 #m2
A2=23992000.0 #m2

alfa=1.0/(A1+A2)

#time interval, seconds, s
DT=86400
#Water density, kg/m3
DH2O=1000.0
#Water specific Heat j/kg C
Cp=4186

s=DT/(2*DZ**2.0)

#Initializing coefficients np.arrays
a=np.zeros(Nlayers-1)
b=np.zeros(Nlayers)
c=np.zeros(Nlayers-1)
d=np.zeros(Nlayers)

for i in xrange(1,N+1):  # 1,2,3,4,5,6,7,8,9,10
#Surface
#b0, c0, d0
b[0]=1+alfa*s*(K+2.0*K+K)
c[0]=-alfa*s*(K+2.0*K+K)
d[0]=((1.0-alfa*s*(K+2.0*K+K))*TempLake[i-1][0])+(alfa*s*(K+2.0*K+K)*TempLake[i-1][1])+alfa*s*(K+K)*((2*Q*DZ)/(Cp*DH2O*K))+alfa*s*(K+K)*((2.0*Q*DZ)/(DH2O*Cp*K)) # 0,1,2,3,4,5,6,7,8,9

#Bottom
#an, bn, dn
a[-1]=-alfa*s*(K+2.0*K+K)
b[-1]=1+alfa*s*(K+2.0*K+K)
d[-1]=alfa*s*TempLake[i-1][-2]*(K+2.0*K+K)+(1-alfa*s*(K+2*K+K))*TempLake[i-1][-1] # 0,1,2,3,4,5,6,7,8,9

#ai, bi, ci, di
a[0]=-alfa*s*(K+K)
a[1]=-alfa*s*(K+K)
a[2]=-alfa*s*(K+K)
a[3]=-alfa*s*(K+K)

b[1]=1.0+alfa*s*(K+2.0*K+K)
b[2]=1.0+alfa*s*(K+2.0*K+K)
b[3]=1.0+alfa*s*(K+2.0*K+K)
b[4]=1.0+alfa*s*(K+2.0*K+K)

c[1]=-alfa*s*(K+K)
c[2]=-alfa*s*(K+K)
c[3]=-alfa*s*(K+K)
c[4]=-alfa*s*(K+K)

d[1]=alfa*s*(K+K)*TempLake[i-1][0]+alfa*s*(K+K)*TempLake[i-1][2]+(1-alfa*s*(K+2.0*K+K))*TempLake[i-1][1]
d[2]=alfa*s*(K+K)*TempLake[i-1][1]+alfa*s*(K+K)*TempLake[i-1][3]+(1-alfa*s*(K+2.0*K+K))*TempLake[i-1][2]
d[3]=alfa*s*(K+K)*TempLake[i-1][2]+alfa*s*(K+K)*TempLake[i-1][4]+(1-alfa*s*(K+2.0*K+K))*TempLake[i-1][3]
d[4]=alfa*s*(K+K)*TempLake[i-1][3]+alfa*s*(K+K)*TempLake[i-1][5]+(1-alfa*s*(K+2.0*K+K))*TempLake[i-1][4]

TempTDMA=TDMAsolver(a, b, c, d)

TempLake[i]=TempTDMA

print TempLake

Z=TempLake

X,Y=np.meshgrid(range(Z.shape[0]+1),range(Z.shape[1]+1))

im = plt.pcolormesh(X,Y,Z.transpose(), cmap='jet')
ax = plt.gca()
ax.set_ylim(ax.get_ylim()[::-1])
plt.colorbar(im, orientation='horizontal')

plt.show()


Results:

TEST 1 : simulation code included above

TEST 2 : simulation equal to Test 1 but lake area is included in km2 instead of m2