fluid dynamics - Buoyancy and natural convection in a stratified water tank

I'm working on a simple calculation-model for a hot-water-storage. At the moment I'm stuck on implementing natural convection.

To calculate the temperature I've implemented a simple 1D-model in python. My tank consists of vertically stacked volumes/cells. If the temperature of one cell is higher than that of the cell above, as shown for cell 5 with a temperature of 50°C, I want to calculate the resulting massflow.

So far I tried to implement the buoyancy-model used in modelica (which is calculating a heat-flow instead of a massflow. That's quite ok...). But I can't get it to work properly. The resulting increase in temperature is just massive (and the units don't fit)... The code for the model can be found here: Modelica buoyancy model
A short example for my calculations (default values are for the mean temperature of T4 and T5: 40°C): $$ \tau = 60\,s, \ \ V = 0.1\,m^3, \ \ T_4 = 30\,^\circ C, \ \ T_5 = 50\,^\circ C, \ \ \rho_{def} = 992.21\,\frac{kg}{m^3}, \ \ c_{p_{def}} = 4178.63\,\frac{J}{kgK} $$ $\Delta T = 20\,K$
$k = \frac{V \rho_{def} c_{p_{def}}}{\tau} = 6910\,\frac{W}{K}$
$Q_{flow} = \dot{Q} = k * \Delta T^2 = 2.7641e6\,\frac{JK}{s}$

$Q = Q_{flow} * \tau = 1.6584e8\,JK$
$\delta T = \frac{Q}{c_{p_{def}}V\rho_{def}} = 400\,K^2$
$sqrt(\delta T) = 20\,K = \Delta T$
So basically the result is my input... Is there a mistake in translating the modelica-code to equations? Or how does modelica get useful results out of this?

My second approach to solve the buoyancy-problem is the use of equations for thermals as in this link at page 171: Plumes and Thermals
But here I'm stuck at calculating the correct buoyancy of cell 5. (I know this sounds simple, but after about 4 days of trying to solve this problem I guess my brain just turned to butter).
When calculation the buoyancy like shown in this image:

And transferring this to my storage-model with cells:
How do I calculate the forces on cell 5? Do I have to consider the density of cells 1 to 3 to calculate the force on the top of cell 5 and cell 6 for the force on the bottom of cell 5?
Or if working with displaced fluid and using the buoyancy forumla $B = \frac{\rho_\infty - \rho_{Cell 5}}{\rho_\infty} g V$: Which density should I use for the surrounding fluid? Am I calculating the density by 50% density of cell 4 and 50% density of cell 6?

Thanks so far! And if anyone got a better idea on how to work on the buoyancy calculation: I'm grateful for your advices!