diff --git a/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py b/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py
index 5d53c01..dc21194 100644
--- a/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py
+++ b/Python/linear_stability/boussinesq/shell/rtc/trace_marginal.py
@@ -130,6 +130,13 @@
 #bc_vel = 0; bc_temp = 0; heating = 0; rratio = 0.6; Pr = 1; rescale = False
 #Ek = 1e-4; Rac = 27.0031; mc = 18
 
+# NS/NS, FT/FT, heating=2, ri/ro=0.35
+#bc_vel = 0; bc_temp = 0
+#heating = 2; alpha=0; beta=1
+#rratio = 0.35; Pr = 0.1; rescale = False
+#Ek = 0.5e-3;  Rac = 18.747; mc = 3
+#res = [32, 32, 0]
+
 # Convert Ekman to Taylor number
 Ta = Ek**-2
 
@@ -149,7 +156,7 @@
     Ek = (Ta*(1.0-rratio)**4)**-0.5
 else:
     Ek = Ta**-0.5
-eq_params = {'ekman':Ek, 'prandtl':Pr, 'rayleigh':Ra, 'r_ratio':rratio, 'heating':heating}
+eq_params = {'ekman':Ek, 'prandtl':Pr, 'rayleigh':Ra, 'r_ratio':rratio, 'heating':heating, 'beta':beta, 'alpha':alpha }
 eq_params.update(model.automatic_parameters(eq_params))
 bcs = {'bcType':model.SOLVER_HAS_BC, 'velocity':bc_vel, 'temperature':bc_temp}
 
@@ -179,7 +186,7 @@ def wave(m):
 marginal_options['save_physical'] = True
 marginal_options['impose_symmetry'] = False
 marginal_options['use_spherical_evp'] = False
-marginal_options['curve_points'] = np.arange(max(1,m-3), m+10, 1)
+marginal_options['curve_points'] = np.arange(max(1,m-3), m+4, 1)
 
 # Compute
 MarginalCurve.compute(gevp_opts, marginal_options)
diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py
index 5d37455..886ee05 100644
--- a/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py
+++ b/Python/quicc_solver/model/boussinesq/shell/rtc/explicit/physical_model.py
@@ -23,7 +23,14 @@ def periodicity(self):
     def nondimensional_parameters(self):
         """Get the list of nondimensional parameters"""
 
-        return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating"]
+        return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating","alpha","beta"]
+
+        # alpha = 1 yields a linear gravity in r, whereas alpha !=1 leads to a alpha*r+(1-alpha)/r^2 type gravity profile
+        # beta = 1 is for pure internal heating, beta =0 is for pure differential heating
+        # heating=0: internal heating;
+        # heating=1: differential heating;
+        # heating=2: differential heating + internal heating (proportional to r)
+        # heating=3: differential heating + internal heating (general form)
 
     def automatic_parameters(self, eq_params):
         """Extend parameters with automatically computable values"""
@@ -236,13 +243,28 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri
         mat = None
         bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_col)
         if field_row == ("velocity","pol") and field_col == ("temperature",""):
-            mat = geo.i4r4(res[0], ri, ro, bc, Ra_eff)
+            alpha = eq_params['alpha']
+            c1 = Ra_eff * alpha
+            c2 = Ra_eff *(1-alpha) * ro**3
+            #mat = geo.i4r4(res[0], ri, ro, bc, Ra_eff)
+            mat = geo.i4r4(res[0], ri, ro, bc, c1) + geo.i4r1(res[0], ri, ro, bc, c2)
 
         elif field_row == ("temperature","") and field_col == ("velocity","pol"):
             if eq_params["heating"] == 0:
                 mat = geo.i2r2(res[0], ri, ro, bc, -bg_eff*l*(l+1.0))
-            else:
+            elif eq_params["heating"] == 1:
                 mat = geo.i2(res[0], ri, ro, bc, -bg_eff*l*(l+1.0))
+            elif eq_params["heating"] == 2 or eq_params["heating"] == 3:
+                beta = eq_params['beta']
+                # assumes length-scale is the depth of the shell
+                c1 = beta*bg_eff/ro # internal heating contribution
+                c2 = ro**2*(1-beta*bg_eff) # differential heating contribution
+                if beta==1/bg_eff: # similar to heating=0
+                    mat = geo.i2r2(res[0], ri, ro, bc, -c1*l*(l+1.0))
+                elif beta==0: # similar to heating=1
+                    mat = geo.i2(res[0], ri, ro, bc, -c2*l*(l+1.0))
+                else:
+                    mat = geo.i2r3(res[0], ri, ro, bc, -c1*l*(l+1.0)) + geo.i2(res[0], ri, ro, bc, -c2*l*(l+1.0))
 
         if mat is None:
             raise RuntimeError("Equations are not setup properly!")
@@ -254,13 +276,21 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr
 
         ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d'])
 
+        Ra_eff, bg_eff = self.nondimensional_factors(eq_params)
+
         mat = None
         bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_col)
         if field_row == ("temperature","") and field_col == field_row:
             if eq_params["heating"] == 0:
                 mat = geo.i2r2(res[0], ri, ro, bc)
-            else:
+            elif eq_params["heating"] == 1:
                 mat = geo.i2r3(res[0], ri, ro, bc)
+            elif eq_params["heating"] == 2 or eq_params["heating"] == 3:
+                beta = eq_params['beta']
+                if beta==1/bg_eff: # same as for heating=0
+                    mat = geo.i2r2(res[0], ri, ro, bc)
+                else:
+                    mat = geo.i2r3(res[0], ri, ro, bc)
 
         if mat is None:
             raise RuntimeError("Equations are not setup properly!")
@@ -277,6 +307,8 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri
 
         ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d'])
 
+        Ra_eff, bg_eff = self.nondimensional_factors(eq_params)
+
         mat = None
         bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_col)
         if field_row == ("velocity","tor") and field_col == field_row:
@@ -288,8 +320,14 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri
         elif field_row == ("temperature","") and field_col == field_row:
             if eq_params["heating"] == 0:
                 mat = geo.i2r2lapl(res[0], ri, ro, l, bc, 1.0/Pr)
-            else:
+            elif eq_params["heating"] == 1:
                 mat = geo.i2r3lapl(res[0], ri, ro, l, bc, 1.0/Pr)
+            elif eq_params["heating"] == 2 or eq_params["heating"] == 3:
+                    beta = eq_params['beta']
+                    if beta==1/bg_eff: # same as for heating=0
+                        mat = geo.i2r2lapl(res[0], ri, ro, l, bc, 1.0/Pr)
+                    else:
+                         mat = geo.i2r3lapl(res[0], ri, ro, l, bc, 1.0/Pr)
 
         if mat is None:
             raise RuntimeError("Equations are not setup properly!")
@@ -304,6 +342,8 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None):
 
         ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d'])
 
+        Ra_eff, bg_eff = self.nondimensional_factors(eq_params)
+        
         mat = None
         bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_row)
         if field_row == ("velocity","tor"):
@@ -315,8 +355,14 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None):
         elif field_row == ("temperature",""):
             if eq_params["heating"] == 0:
                 mat = geo.i2r2(res[0], ri, ro, bc)
-            else:
+            elif eq_params["heating"] == 1:
                 mat = geo.i2r3(res[0], ri, ro, bc)
+            elif eq_params["heating"] == 2 or eq_params["heating"] == 3:
+                beta = eq_params['beta']
+                if beta == 1/bg_eff: # same as for heating=0
+                    mat = geo.i2r2(res[0], ri, ro, bc)
+                else:
+                    mat = geo.i2r3(res[0], ri, ro, bc)
 
         if mat is None:
             raise RuntimeError("Equations are not setup properly!")
@@ -361,5 +407,11 @@ def nondimensional_factors(self, eq_params):
             # (R_o - R_i) rescaling
             Ra_eff = (Ra*T/ro)
             bg_eff = ro**2*rratio
+        elif eq_params['heating'] == 2:
+            Ra_eff = (Ra*T/ro)
+            bg_eff = 1
+        elif eq_params['heating'] == 3:
+            Ra_eff = (Ra*T/ro)
+            bg_eff = 1/(1-rratio**3)
 
         return (Ra_eff, bg_eff)
diff --git a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py
index 9ff48d6..d2654c0 100644
--- a/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py
+++ b/Python/quicc_solver/model/boussinesq/shell/rtc/implicit/physical_model.py
@@ -22,7 +22,14 @@ def periodicity(self):
     def nondimensional_parameters(self):
         """Get the list of nondimensional parameters"""
 
-        return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating"]
+        return ["ekman", "prandtl", "rayleigh", "r_ratio", "heating","alpha","beta"]
+
+        # alpha = 1 yields a linear gravity in r, whereas alpha !=1 leads to a alpha*r+(1-alpha)/r^2 type gravity profile
+        # beta = 1 is for pure internal heating, beta=0 is for pure differential heating
+        # heating=0: internal heating;   
+        # heating=1: differential heating; 
+        # heating=2: differential heating + internal heating (proportional to r)
+        # heating=3: differential heating + internal heating (general form)
 
     def automatic_parameters(self, eq_params):
         """Extend parameters with automatically computable values"""
@@ -192,12 +199,16 @@ def convert_bc(self, eq_params, eigs, bcs, field_row, field_col):
                         bc = {0:-22, 'rt':0}
                     elif field_col == ("velocity","pol"):
                         bc = {0:-41, 'rt':0}
+                    elif field_col == ("temperature",""):
+                        bc = {0:-21, 'rt':0}
 
                 else:
                     if field_row == ("velocity","tor") and field_col == ("velocity","tor"):
                         bc = {0:22}
                     elif field_row == ("velocity","pol") and field_col == ("velocity","pol"):
                         bc = {0:41}
+                    elif field_row == ("temperature","") and field_col == ("temperature",""):
+                            bc = {0:21}
 
             # Set LHS galerkin restriction
             if self.use_galerkin:
@@ -258,8 +269,19 @@ def explicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri
         if field_row == ("temperature","") and field_col == ("velocity","pol"):
             if eq_params["heating"] == 0:
                 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, -bg_eff, with_sh_coeff = 'laplh', restriction = restriction)
-            else:
+            elif eq_params["heating"] == 1:
                 mat = geo.i2(res[0], ri, ro, res[1], m, bc, -bg_eff, with_sh_coeff = 'laplh', restriction = restriction)
+            elif eq_params["heating"] == 2 or eq_params["heating"] == 3:
+                beta = eq_params['beta']
+                # assumes length-scale is the depth of the shell
+                c1 = -beta*bg_eff/ro # internal heating contribution
+                c2 = -ro**2*(1-beta*bg_eff) # differential heating contribution
+                if beta==1/bg_eff: # similar to heating=0
+                    mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction)
+                elif beta==0: # similar to heating==1
+                    mat = geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction)
+                else:
+                    mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction)
 
         if mat is None:
             raise RuntimeError("Equations are not setup properly!")
@@ -275,13 +297,21 @@ def nonlinear_block(self, res, eq_params, eigs, bcs, field_row, field_col, restr
 
         ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d'])
 
+        Ra_eff, bg_eff = self.nondimensional_factors(eq_params)
+
         mat = None
         bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_col)
         if field_row == ("temperature","") and field_col == field_row:
             if eq_params["heating"] == 0:
                 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction)
-            else:
+            elif eq_params["heating"] == 1:
                 mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction)
+            elif eq_params["heating"] == 2 or eq_params["heating"] == 3:
+                beta = eq_params['beta']
+                if beta==1/bg_eff: # similar to heating=0
+                    mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction)
+                else:
+                    mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction)
 
         if mat is None:
             raise RuntimeError("Equations are not setup properly!")
@@ -331,7 +361,11 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri
                 mat = mat + geo.i4r4lapl(res[0], ri, ro, res[1], m, bc, 1j*m*T/c_dt, with_sh_coeff = 'invlaplh', l_zero_fix = 'zero', restriction = restriction)
 
             elif field_col == ("temperature",""):
-                mat = geo.i4r4(res[0], ri, ro, res[1], m, bc, -Ra_eff/c_dt**2, l_zero_fix = 'zero', restriction = restriction)
+                alpha = eq_params['alpha']
+                c1 = -Ra_eff * alpha /c_dt**2
+                c2 = -Ra_eff *(1-alpha) * ro**3 /c_dt**2
+                #mat = geo.i4r4(res[0], ri, ro, res[1], m, bc, -Ra_eff/c_dt**2, l_zero_fix = 'zero', restriction = restriction)
+                mat = geo.i4r4(res[0], ri, ro, res[1], m, bc, c1, l_zero_fix = 'zero', restriction = restriction) + geo.i4r1(res[0], ri, ro, res[1], m, bc, c2, l_zero_fix = 'zero', restriction = restriction)
 
         elif field_row == ("temperature",""):
             if field_col == ("velocity","tor"):
@@ -340,9 +374,20 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri
             elif field_col == ("velocity","pol"):
                 if self.linearize:
                     if eq_params["heating"] == 0:
-                        mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, bg_eff/Pr, with_sh_coeff = 'laplh', restriction = restriction)
-                    else:
-                        mat = geo.i2(res[0], ri, ro, res[1], m, bc, bg_eff/Pr, with_sh_coeff = 'laplh', restriction = restriction)
+                        mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, bg_eff, with_sh_coeff = 'laplh', restriction = restriction)
+                    elif eq_params["heating"] == 1:
+                        mat = geo.i2(res[0], ri, ro, res[1], m, bc, bg_eff, with_sh_coeff = 'laplh', restriction = restriction)
+                    elif eq_params["heating"] == 2 or eq_params["heating"] == 3:
+                        beta = eq_params['beta']
+                        # assumes length-scale is the depth of the shell
+                        c1 = beta*bg_eff/ro # internal heating contribution
+                        c2 = ro**2*(1-beta*bg_eff) # differential heating contribution
+                        if beta==1/bg_eff: # similar to heating=0
+                            mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction)
+                        elif beta==0:  # similar to heating=1
+                            mat = geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction)
+                        else:
+                            mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, c1, with_sh_coeff = 'laplh', restriction = restriction) + geo.i2(res[0], ri, ro, res[1], m, bc, c2, with_sh_coeff = 'laplh', restriction = restriction)
 
                 else:
                     mat = geo.zblk(res[0], ri, ro, res[1], m, bc)
@@ -350,8 +395,14 @@ def implicit_block(self, res, eq_params, eigs, bcs, field_row, field_col, restri
             elif field_col == ("temperature",""):
                 if eq_params["heating"] == 0:
                     mat = geo.i2r2lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction)
-                else:
+                elif eq_params["heating"] == 1:
                     mat = geo.i2r3lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction)
+                elif eq_params["heating"] == 2 or eq_params["heating"] == 3:
+                    beta = eq_params['beta']
+                    if beta==1/bg_eff: # same as for heating=0
+                        mat = geo.i2r2lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction)
+                    else:
+                        mat = geo.i2r3lapl(res[0], ri, ro, res[1], m, bc, 1.0/(Pr*c_dt), restriction = restriction)
 
         if mat is None:
             raise RuntimeError("Equations are not setup properly!")
@@ -367,6 +418,8 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None):
 
         ri, ro = (self.automatic_parameters(eq_params)['lower1d'], self.automatic_parameters(eq_params)['upper1d'])
 
+        Ra_eff, bg_eff = self.nondimensional_factors(eq_params)
+
         mat = None
         bc = self.convert_bc(eq_params,eigs,bcs,field_row,field_row)
         if field_row == ("velocity","tor"):
@@ -378,8 +431,14 @@ def time_block(self, res, eq_params, eigs, bcs, field_row, restriction = None):
         elif field_row == ("temperature",""):
             if eq_params["heating"] == 0:
                 mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction)
-            else:
+            elif eq_params["heating"] == 1:
                 mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction)
+            elif eq_params["heating"] == 2 or eq_params["heating"] == 3:
+                beta = eq_params['beta']
+                if beta == 1/bg_eff: # same as for heating=0
+                    mat = geo.i2r2(res[0], ri, ro, res[1], m, bc, restriction = restriction)
+                else:
+                    mat = geo.i2r3(res[0], ri, ro, res[1], m, bc, restriction = restriction)
 
         if mat is None:
             raise RuntimeError("Equations are not setup properly!")
@@ -440,9 +499,14 @@ def nondimensional_factors(self, eq_params):
         elif eq_params['heating'] == 1:
             # (R_o - R_i) rescaling
             Ra_eff = (Ra*T/ro)
-            bg_eff = ro**2*rratio
-
+            bg_eff = ro**2*rratio       
         # Rescale Rayleigh by E^{-4/3}
         # Ra_eff = Ra_eff*T**(1./3.)
+        elif eq_params['heating'] == 2:
+            Ra_eff = (Ra*T/ro)
+            bg_eff = 1
+        elif eq_params['heating'] == 3:
+            Ra_eff = (Ra*T/ro)
+            bg_eff = 1/(1-rratio**3)
 
         return (Ra_eff, bg_eff)
diff --git a/Readme.md b/Readme.md
index f8144bf..91a942d 100644
--- a/Readme.md
+++ b/Readme.md
@@ -1,15 +1,83 @@
 # Boussinesq rotating thermal convection(RTC) in a spherical shell
 
+
 The model equations are
 
-![Model equations](https://latex.codecogs.com/svg.latex?%5Cinline%20%5Cbegin%7Balign*%7D%20%5Cleft%28%5Cpartial_t%20-%20%5CDelta%5Cright%29%5Cmathbf%7Bu%7D%20%26%20%3D%20%5Cmathbf%7Bu%7D%20%5Ctimes%20%5Cleft%28%20%5Cnabla%20%5Ctimes%20%5Cmathbf%7Bu%7D%5Cright%29%20+%20%5Cfrac%7BRa%7D%7BPr%7D%20%5CTheta%20%5Cmathbf%7Br%7D%20-%20%5Cnabla%5CPi%20%5C%5C%5B0.6cm%5D%20%5Cleft%28%5Cpartial_t%20-%20%5Cfrac%7B1%7D%7BPr%7D%5CDelta%5Cright%29%5CTheta%20%26%20%3D%20%5Cfrac%7BS%7D%7BPr%7D%20-%20%5Cmathbf%7Bu%7D%5Ccdot%5Cnabla%5CTheta%5C%5C%5B0.6cm%5D%20%5Cnabla%5Ccdot%5Cmathbf%7Bu%7D%20%26%20%3D%200%20%5Cend%7Balign*%7D)
+$$
+(\partial_t -\Delta){\bf u}= {\bf u}\times(\nabla\times{\bf u})- \frac{1}{E}\hat{\bf z}\times{\bf u}  -\nabla \Pi + \frac{Ra}{E} \left(\frac{d}{r_o}\right)\Theta{\bf r}\\
+$$
+
+$$
+(\partial_t -\frac{1}{Pr}\Delta)\Theta = -{\bf u}\cdot\nabla\Theta + \mathcal{H}{\bf u}\cdot{\bf r}
+$$
 
 with the parameters defined as
 
-![Nondimensional parameters](https://latex.codecogs.com/svg.latex?%5Cinline%20%5Cbegin%7Balign*%7D%20Pr%20%26%20%3D%20%5Cfrac%7B%5Cnu%7D%7B%5Ckappa%7D%5C%5C%20Ra%20%26%20%3D%20%5Cfrac%7Bg%20%5Calpha%20%5Cbeta%20r_o%5E4%7D%7B%5Cnu%5Ckappa%7D%20%5Cend%7Balign*%7D)
+$$
+Pr = \frac{\nu}{\kappa};\quad E=\frac{\nu}{2\Omega d^2};\quad Ra=\frac{ \alpha g_o \mathcal{T} d}{2\nu\Omega}
+$$
+
+where $\nu$, $\kappa$, $\alpha$ are the kinematic viscosity, thermal conductivity, and thermal expansion coefficents; $\Omega$ is the rotation rate; $r_o$ is the outer radius of the spherical shell and $d$ is its thickness; $g_o$ is defined via the gravitational acceleration, ${\bf g} = -g_o {\bf r}/r_o$. $\mathcal{T}$ is a characteristic temperature and $\mathcal{H}$ is the effective temperature background. The spherical shell has an inner radius of $r_i$, outeer radius $r_o$ and thickness $d=r_o-r_i$.
+
+The model equations are non-dimensionalised using $d$ as the characteristic length-scale and the viscous time-scale $\nu/d^2$. Note that in the model equation, the factor $d/r_o$ as written in the buoyancy term, involves dimensional quantities. If the outer radius $r_o$ is already non-dimensional, that factor is simply $1/r_o$.
+
+The definition of both $\mathcal{T}$ and $\mathcal{H}$ depend on the mode of thermal driving (differential or internal heating, controlled by the parameter `heating`) and on the boundary conditions (fixed-temperature or fixed-flux):
+
+- For fixed temperature boundary conditions, the background temperature is as given in Eq. 1.6 of Dormy et al., 2004 (DOI: 10.1017/S0022112003007316). Then $\mathcal{T} = \Delta T = T_i - T_o$ (where $T$ is the temperature of the resulting reference state) and:
+    + $\mathcal{H}=$`bg_eff` for uniform internal heating (`heating=0`)
+    + $\mathcal{H}=$ `bg_eff` $/r^3$ for differential heating (`heating=1`)  
+where `bg_eff` is a pre-factor, function of $r_o$ and $d$.
+
+- **The fixed-flux case is not currently implemented cleanly**. For fixed-flux at either one of the boundaries, the relevant thermal scale is defined by a thermal gradient $\beta$, so that the characteristic temperature is $\mathcal{T}=\beta d$. The exact definition of $\beta$ (and therefore of the Rayleigh number) depends on the mode of heating, but still:
+    + $\mathcal{H}=$`bg_eff` for uniform internal heating (`heating=0`)
+    + $\mathcal{H}=$ `bg_eff` $/r^3$ for differential heating (`heating=1`)  
+**Currently the correct** `bg_eff` **is not implemented for fixed-flux conditions, but the code can still be used to calculate the onset of convection. The resulting Rayleigh number needs, however, to be corrected to account for this discrepancy.**
+
+
+**At present, only uniform heating is implemented in the internally heated case.**
+
+The system is solved using a Toroidal/Poloidal decomposition of the velocity
+
+$$
+{\bf u} = \nabla\times T {\bf r} + \nabla\times\nabla\times P {\bf r} .
+$$
+
+
+## General background profiles
+
+**The following functionalities are currently implemented only for the Explicit models, and not for the SplitEqaution formalism.**
+
+The parameters $\alpha$ (`alpha`) and $\beta$ (`beta`) can be used to set up general gravity and background temperature profiles. The following are currently implemented:
+
+### Gravity:
+The following general gravity profile is implemented:
+
+$$
+{\bf g} = -g_o\left[\left(\frac{r}{r_o}\right)\alpha + \left(\frac{r_o}{r}\right)^2(1-\alpha)\right]\hat{\bf r}
+$$
+
+where, with $\alpha=1$ a linear profile is obtained, adequate for an inner core of the same density as the outer core.
+
+**This is the feature not implemented in the SplitEqaution formalism.**
+
+### Temperature:
+
+General background temperature profiles, $T$, containing a mix of both internal and differential heating are implemented.
+
+There are two options, chosed via the parameter `heating`:
+
+- `heating`=2 implements a purely linear internal heating component, resulting in:
+
+$$
+T = T'_o\left[\left(\frac{r}{r_o}\right)\beta + \left(\frac{r_o}{r}\right)^2(1-\beta)\right]
+$$
+
+- `heating`=3 implements a fully general form for the  internal heating component, resulting in:
 
-The system is solved using a Toroidal/Poloidal decomposition of the velocity ![u](https://latex.codecogs.com/svg.latex?%5Cinline%20%5Cmathbf%7Bu%7D):
+$$
+\frac{dT}{dr} = T'_o\left[\frac{1}{1-\chi^3}\left(\frac{r}{r_o}\right)\beta + \left(\frac{r_o}{r}\right)^2\left(1-\frac{\beta}{1-\chi^3}\right)\right]
+$$
 
-![Toroidal/Poloidal decomposition](https://latex.codecogs.com/svg.latex?%5Cinline%20%5Cmathbf%7Bu%7D%3D%5Cmathbf%7B%5Cnabla%7D%5Ctimes%20T%20%5Cmathbf%7Br%7D%20+%20%5Cmathbf%7B%5Cnabla%7D%5Ctimes%5Cmathbf%7B%5Cnabla%7D%5Ctimes%20P%20%5Cmathbf%7Br%7D)
+In both formulations, `heating`=1 is the pure internal heating case, and `heating`=0 is the differential heating case.
 
-Note: To access the equations in codecogs online LaTeX editor replace "https://latex.codecogs.com/svg.latex?" with "https://www.codecogs.com/eqnedit.php?latex=".
+In both cases, tempreature is nondimensionalised via $T'_o$, the temperature gradient at the outer boundary, $r_o$. In the above definition of the Rayleigh number, therefore: $\mathcal{T} = D T'_o$.
\ No newline at end of file
diff --git a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp
index 4541c18..ed01894 100644
--- a/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp
+++ b/src/Model/Boussinesq/Shell/RTC/Explicit/ModelBackend.cpp
@@ -33,8 +33,12 @@
 #include "QuICC/PhysicalNames/Temperature.hpp"
 #include "QuICC/NonDimensional/Prandtl.hpp"
 #include "QuICC/NonDimensional/Rayleigh.hpp"
+#include "QuICC/NonDimensional/Alpha.hpp"
+#include "QuICC/NonDimensional/Beta.hpp"
 #include "QuICC/NonDimensional/Ekman.hpp"
 #include "QuICC/NonDimensional/Heating.hpp"
+#include "QuICC/NonDimensional/Alpha.hpp"
+#include "QuICC/NonDimensional/Beta.hpp"
 #include "QuICC/NonDimensional/RRatio.hpp"
 #include "QuICC/NonDimensional/Lower1d.hpp"
 #include "QuICC/NonDimensional/Upper1d.hpp"
@@ -49,6 +53,7 @@
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y2SphLapl.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3SphLapl.hpp"
+#include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y1.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl2.hpp"
@@ -231,17 +236,32 @@ namespace Explicit {
       {
          auto Pr = nds.find(NonDimensional::Prandtl::id())->second->value();
          auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value();
+         auto beta = nds.find(NonDimensional::Beta::id())->second->value();
+         auto bg = effectiveBg(nds);
 
          if(heatingMode == 0)
          {
             SparseSM::Chebyshev::LinearMap::I2Y2SphLapl spasm(nN, nN, ri, ro, l);
             decMat.real() = (1.0/Pr)*spasm.mat();
          }
-         else
+         else if(heatingMode == 1)
          {
             SparseSM::Chebyshev::LinearMap::I2Y3SphLapl spasm(nN, nN, ri, ro, l);
             decMat.real() = (1.0/Pr)*spasm.mat();
          }
+         else if(heatingMode == 2 || heatingMode == 3)
+         {
+            if(beta == 1/bg)
+            {
+               SparseSM::Chebyshev::LinearMap::I2Y2SphLapl spasm(nN, nN, ri, ro, l);
+               decMat.real() = (1.0/Pr)*spasm.mat();
+            }
+            else
+            {
+               SparseSM::Chebyshev::LinearMap::I2Y3SphLapl spasm(nN, nN, ri, ro, l);
+               decMat.real() = (1.0/Pr)*spasm.mat();
+            }
+         }
       }
       else
       {
@@ -280,17 +300,32 @@ namespace Explicit {
       else if(fieldId == std::make_pair(PhysicalNames::Temperature::id(), FieldComponents::Spectral::SCALAR))
       {
          auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value();
+         auto beta = nds.find(NonDimensional::Beta::id())->second->value();
+         auto bg = effectiveBg(nds);
 
          if(heatingMode == 0)
          {
             SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro);
             decMat.real() = spasm.mat();
          }
-         else
+         else if(heatingMode == 1)
          {
             SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro);
             decMat.real() = spasm.mat();
          }
+         else if(heatingMode == 2 || heatingMode == 3)
+         {
+            if(beta == 1/bg)
+            {
+               SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro);
+               decMat.real() = spasm.mat();
+            }
+            else
+            {
+               SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro);
+               decMat.real() = spasm.mat();
+            }
+         }
       }
    }
 
@@ -423,6 +458,8 @@ namespace Explicit {
       auto ri = nds.find(NonDimensional::Lower1d::id())->second->value();
       auto ro = nds.find(NonDimensional::Upper1d::id())->second->value();
 
+      auto alpha = nds.find(NonDimensional::Alpha::id())->second->value();
+
       int s = this->nBc(rowId);
       // Explicit linear operator
       if(opId == ModelOperator::ExplicitLinear::id())
@@ -431,22 +468,45 @@ namespace Explicit {
                colId == std::make_pair(PhysicalNames::Temperature::id(),FieldComponents::Spectral::SCALAR))
          {
             auto Ra = effectiveRa(nds);
+            MHDFloat c1, c2;
+            c1 = Ra*alpha;
+            c2 = Ra*(1.-alpha)*(ro*ro*ro);
 
             if(this->useSplitEquation())
             {
-               SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro);
-               decMat.real() = Ra*spasm.mat();
+               if(alpha == 1) // linear gravity
+               {
+                  SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro);
+                  decMat.real() = Ra*spasm.mat();
+               }
+               else //general gravity
+               {
+                  // To fully implement the splitEquation version , we need to update the nonlinear operator
+                  // work in progress
+                  SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro); 
+                  decMat.real() = Ra*spasm.mat();                  
+               }
             }
             else
             {
-               SparseSM::Chebyshev::LinearMap::I4Y4 spasm(nN, nN, ri, ro);
-               decMat.real() = Ra*spasm.mat();
+               if(alpha == 1) //linear gravity
+               {
+                  SparseSM::Chebyshev::LinearMap::I4Y4 spasm(nN, nN, ri, ro);
+                  decMat.real() = Ra*spasm.mat();
+               }
+               else // general gravity
+               {
+                  SparseSM::Chebyshev::LinearMap::I4Y4 i4r4(nN, nN, ri, ro);
+                  SparseSM::Chebyshev::LinearMap::I4Y1 i4r1(nN, nN, ri, ro);
+                  decMat.real() = c1*i4r4.mat() + c2*i4r1.mat();
+               }
             }
          }
          else if(rowId == std::make_pair(PhysicalNames::Temperature::id(), FieldComponents::Spectral::SCALAR) && 
                colId == std::make_pair(PhysicalNames::Velocity::id(),FieldComponents::Spectral::POL))
          {
             auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value();
+            auto beta = nds.find(NonDimensional::Beta::id())->second->value();
             auto bg = effectiveBg(nds);
             auto dl = static_cast<MHDFloat>(l);
             auto ll1 = dl*(dl + 1.0);
@@ -456,11 +516,33 @@ namespace Explicit {
                SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro);
                decMat.real() = -bg*ll1*spasm.mat();
             }
-            else
+            else if(heatingMode == 1)
             {
                SparseSM::Chebyshev::LinearMap::I2 spasm(nN, nN, ri, ro);
                decMat.real() = -bg*ll1*spasm.mat();
             }
+            else if(heatingMode == 2 || heatingMode == 3)
+            {
+               MHDFloat c1, c2;
+               c1 = beta*bg/ro;
+               c2 = ro*ro*(1-beta*bg);
+               if(beta == 1/bg) 
+               {
+                  SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro);
+                  decMat.real() = -c1*ll1*spasm.mat();
+               }
+               else if(beta == 0)
+               {
+                  SparseSM::Chebyshev::LinearMap::I2 spasm(nN, nN, ri, ro);
+                  decMat.real() = -c2*ll1*spasm.mat();
+               }
+               else
+               {
+                  SparseSM::Chebyshev::LinearMap::I2Y3 i2r3(nN, nN, ri, ro);
+                  SparseSM::Chebyshev::LinearMap::I2 i2(nN, nN, ri, ro);
+                  decMat.real() = -c1*ll1*i2r3.mat() -c2*ll1*i2.mat();
+               }
+            }
          }
          else
          {
@@ -474,17 +556,32 @@ namespace Explicit {
          if(rowId == std::make_pair(PhysicalNames::Temperature::id(), FieldComponents::Spectral::SCALAR) && rowId == colId)
          {
             auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value();
+            auto beta = nds.find(NonDimensional::Beta::id())->second->value();
+            auto bg = effectiveBg(nds);
 
             if(heatingMode == 0)
             {
                SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro);
                decMat.real() = spasm.mat();
             }
-            else
+            else if(heatingMode == 1)
             {
                SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro);
                decMat.real() = spasm.mat();
             }
+            else if(heatingMode == 2 || heatingMode == 3)
+            {
+               if(beta==1/bg)
+               {
+                  SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro);
+                  decMat.real() = spasm.mat();
+               }
+               else
+               {
+                  SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro);
+                  decMat.real() = spasm.mat();
+               }
+            }
          }
          else
          {
diff --git a/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp b/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp
index f8645b5..02fec8b 100644
--- a/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp
+++ b/src/Model/Boussinesq/Shell/RTC/IRTCBackend.cpp
@@ -30,6 +30,8 @@
 #include "QuICC/PhysicalNames/Temperature.hpp"
 #include "QuICC/NonDimensional/Prandtl.hpp"
 #include "QuICC/NonDimensional/Rayleigh.hpp"
+#include "QuICC/NonDimensional/Alpha.hpp"
+#include "QuICC/NonDimensional/Beta.hpp"
 #include "QuICC/NonDimensional/Ekman.hpp"
 #include "QuICC/NonDimensional/Heating.hpp"
 #include "QuICC/NonDimensional/RRatio.hpp"
@@ -46,6 +48,7 @@
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y2SphLapl.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3SphLapl.hpp"
+#include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y1.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl2.hpp"
@@ -91,6 +94,8 @@ namespace RTC {
       std::vector<std::string> names = {
          NonDimensional::Prandtl().tag(),
          NonDimensional::Rayleigh().tag(),
+	 NonDimensional::Alpha().tag(),
+	 NonDimensional::Beta().tag(),
          NonDimensional::Ekman().tag(),
          NonDimensional::Heating().tag(),
          NonDimensional::RRatio().tag()
@@ -151,6 +156,7 @@ namespace RTC {
       auto ro = nds.find(NonDimensional::Upper1d::id())->second->value();
       auto rratio = nds.find(NonDimensional::RRatio::id())->second->value();
       auto heatingMode = nds.find(NonDimensional::Heating::id())->second->value();
+      auto beta = nds.find(NonDimensional::Beta::id())->second->value();
 
       if(ro == 1.0)
       {
@@ -167,6 +173,15 @@ namespace RTC {
       {
          effBg = ro*ro*rratio;
       }
+      // gap with mixed heating
+      else if(heatingMode == 2)
+      {
+         effBg = 1.0;
+      }
+      else if(heatingMode == 3)
+      {
+         effBg = 1.0 / (1-rratio*rratio*rratio);
+      }
 
       return effBg;
    }
diff --git a/src/Model/Boussinesq/Shell/RTC/IRTCModel.cpp b/src/Model/Boussinesq/Shell/RTC/IRTCModel.cpp
index 3385315..fb449da 100644
--- a/src/Model/Boussinesq/Shell/RTC/IRTCModel.cpp
+++ b/src/Model/Boussinesq/Shell/RTC/IRTCModel.cpp
@@ -28,6 +28,8 @@
 #include "QuICC/PhysicalNames/VorticityZ.hpp"
 #include "QuICC/NonDimensional/Ekman.hpp"
 #include "QuICC/NonDimensional/Rayleigh.hpp"
+#include "QuICC/NonDimensional/Alpha.hpp"
+#include "QuICC/NonDimensional/Beta.hpp"
 #include "QuICC/NonDimensional/Prandtl.hpp"
 #include "QuICC/NonDimensional/RRatio.hpp"
 #include "QuICC/NonDimensional/Heating.hpp"
diff --git a/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp b/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp
index 3cecee3..bce8b18 100644
--- a/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp
+++ b/src/Model/Boussinesq/Shell/RTC/Implicit/ModelBackend.cpp
@@ -32,6 +32,8 @@
 #include "QuICC/PhysicalNames/Temperature.hpp"
 #include "QuICC/NonDimensional/Prandtl.hpp"
 #include "QuICC/NonDimensional/Rayleigh.hpp"
+#include "QuICC/NonDimensional/Alpha.hpp"
+#include "QuICC/NonDimensional/Beta.hpp"
 #include "QuICC/NonDimensional/Ekman.hpp"
 #include "QuICC/NonDimensional/Heating.hpp"
 #include "QuICC/NonDimensional/RRatio.hpp"
@@ -52,6 +54,7 @@
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I2Y3SphLapl.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y3.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4D1.hpp"
+#include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y1.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl.hpp"
 #include "QuICC/SparseSM/Chebyshev/LinearMap/I4Y4SphLapl2.hpp"
@@ -488,7 +491,7 @@ namespace Implicit {
                   SparseSM::Chebyshev::LinearMap::I2Y2SphLapl spasm(nN, nN, ri, ro, l);
                   bMat = (1.0/Pr)*spasm.mat();
                }
-               else
+               else if(heatingMode == 1)
                {
                   SparseSM::Chebyshev::LinearMap::I2Y3SphLapl spasm(nN, nN, ri, ro, l);
                   bMat = (1.0/Pr)*spasm.mat();
@@ -605,7 +608,7 @@ namespace Implicit {
                SparseSM::Chebyshev::LinearMap::I2Y2 spasm(nN, nN, ri, ro);
                bMat = spasm.mat();
             }
-            else
+            else if(heatingMode == 1)
             {
                SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro);
                bMat = spasm.mat();
@@ -852,7 +855,7 @@ namespace Implicit {
                   }
                   //this->addBlock(decMat.real(), spasm.mat(), rowShift, colShift, -bg*ll1);
                }
-               else
+               else if(heatingMode == 1)
                {
                   SparseSM::Chebyshev::LinearMap::I2 spasm(nN, nN, ri, ro);
                   bMat = -bg*ll1*spasm.mat();
@@ -904,7 +907,7 @@ namespace Implicit {
                   }
                   //this->addBlock(decMat.real(), spasm.mat(), rowShift, colShift);
                }
-               else
+               else if(heatingMode == 1)
                {
                   SparseSM::Chebyshev::LinearMap::I2Y3 spasm(nN, nN, ri, ro);
                   bMat = spasm.mat();