Two-layer Benard Rayleigh Natural Convection

I have defined the gravity force =3w(i)rho(T-Tref)gbeta/(Th-Tc)
for two layers should I define the force due to the dnesity variaition cuased by the croosing interface region?
::::::::::::::::::::;;;
Also, May I have some good refrences?. I have built my own code and I can share it with any person to help each others.