diff --git a/constants.py b/constants.py
index 45340ec7..2e118946 100644
--- a/constants.py
+++ b/constants.py
@@ -17,6 +17,7 @@
 roughness_method = 'Moelg12'                    # possibilities: 'Moelg12'
 saturation_water_vapour_method = 'Sonntag90'    # possibilities: 'Sonntag90'
 thermal_conductivity_method = 'bulk'		    # possibilities: 'bulk', 'empirical'
+specific_heat_capacity_method = 'bulk'          # possibilities: 'bulk', 'Yen81'
 sfc_temperature_method = 'SLSQP'                # possibilities: 'L-BFGS-B', 'SLSQP'(faster), 'Newton' (Secant, fastest)'
 
 # WRF_X_CSPY: for efficiency and consistency
diff --git a/cosipy/cpkernel/node.py b/cosipy/cpkernel/node.py
index 88a6abce..3bf8576c 100644
--- a/cosipy/cpkernel/node.py
+++ b/cosipy/cpkernel/node.py
@@ -131,7 +131,15 @@ def get_layer_specific_heat(self):
             cp : float
                 Specific heat [:math:`J~kg^{-1}~K^{-1}`]
         """
-        return self.get_layer_ice_fraction()*spec_heat_ice + self.get_layer_air_porosity()*spec_heat_air + self.get_layer_liquid_water_content()*spec_heat_water
+
+        methods_allowed = ['bulk','Yen81']
+        if specific_heat_capacity_method == 'bulk':    
+            cp = self.get_layer_ice_fraction() * spec_heat_ice + self.get_layer_air_porosity() * spec_heat_air + self.get_layer_liquid_water_content() * spec_heat_water
+        elif specific_heat_capacity_method == 'Yen81':    # empirical parameterisation according to Yen, 1981
+            cp = 152.2 + 7.122 * self.get_layer_temperature()
+        else:
+            raise ValueError("Specific heat method = \"{:s}\" is not allowed, must be one of {:s}".format(specific_heat_capacity_method, ", ".join(methods_allowed)))
+        return cp
 
     def get_layer_liquid_water_content(self):
         """ Returns the liquid water content of the node.