diff --git a/src/reaction_modules/asm3_class.js b/src/reaction_modules/asm3_class.js
index ab272dc..d228619 100644
--- a/src/reaction_modules/asm3_class.js
+++ b/src/reaction_modules/asm3_class.js
@@ -69,6 +69,28 @@ class ASM3 {
       i_cNH:    1/14, // charge per S_NH                [mole H+ g-1 NH3-N]
       i_cNO:   -1/14  // charge per S_NO                [mole H+ g-1 NO3-N]
     };
+
+    /**
+     * Temperature theta parameters for ASM3.
+     * These parameters are used to adjust reaction rates based on temperature.
+     * @property {Object} temp_params - Temperature theta parameters
+     */
+    this.temp_params = {
+      // Hydrolysis
+      theta_H:        this._compute_theta(2, 3, 10, 20),
+      // Heterotrophs
+      theta_STO:      this._compute_theta(2.5, 5, 10, 20),
+      theta_mu_H:     this._compute_theta(1, 2, 10, 20),
+      theta_b_H_O:    this._compute_theta(0.1, 0.2, 10, 20),
+      theta_b_H_NO:   this._compute_theta(0.05, 0.1, 10, 20),
+      theta_b_STO_O:  this._compute_theta(0.1, 0.2, 10, 20),
+      theta_b_STO_NO: this._compute_theta(0.05, 0.1, 10, 20),
+      // Autotrophs
+      theta_mu_A:     this._compute_theta(0.35, 1, 10, 20),
+      theta_b_A_O:    this._compute_theta(0.05, 0.15, 10, 20),
+      theta_b_A_NO:   this._compute_theta(0.02, 0.05, 10, 20)
+    };
+
     this.stoi_matrix = this._initialise_stoi_matrix();
   }
 
@@ -103,7 +125,7 @@ class ASM3 {
    * @param {number} K - Half-saturation constant for the reaction species.
    * @returns {number} - Monod equation rate value for the given concentration and half-saturation constant.
    */
-  _monod(c, K){
+  _monod(c, K) {
     return c / (K + c);
   }
 
@@ -113,50 +135,76 @@ class ASM3 {
    * @param {number} K - Half-saturation constant for the reaction species. 
    * @returns {number} - Inverse Monod equation rate value for the given concentration and half-saturation constant.
    */
-  _inv_monod(c, K){
+  _inv_monod(c, K) {
     return K / (K + c);
   }
 
   /**
-   * Computes the reaction rates for each process reaction based on the current state.
+   * Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter.
+   * @param {number} k - Rate constant at 20 degrees Celcius.
+   * @param {number} theta - Theta parameter.
+   * @param {number} T - Temperature in Celcius.
+   * @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation.
+   */
+  _arrhenius(k, theta, T) {
+    return k * Math.exp(theta*(T-20));
+  }
+
+  /**
+   * Computes the temperature theta parameter based on two rate constants and their corresponding temperatures.
+   * @param {number} k1 - Rate constant at temperature T1.
+   * @param {number} k2 - Rate constant at temperature T2.
+   * @param {number} T1 - Temperature T1 in Celcius.
+   * @param {number} T2 - Temperature T2 in Celcius.
+   * @returns {number} - Theta parameter.
+   */
+  _compute_theta(k1, k2, T1, T2) {
+    return Math.log(k1/k2)/(T1-T2);
+  }
+
+  /**
+   * Computes the reaction rates for each process reaction based on the current state and temperature.
    * @param {Array} state - State vector containing concentrations of reaction species.
+   * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
    * @returns {Array} - Reaction rates for each process reaction.
    */
-  compute_rates(state) {
+  compute_rates(state, T = 20) {
     // state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
     const rates = Array(12);
     const [S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS] = state;
     const { k_H, K_X, k_STO, nu_NO, K_O, K_NO, K_S, K_STO, mu_H_max, K_NH, K_HCO, b_H_O, b_H_NO, b_STO_O, b_STO_NO, mu_A_max, K_A_NH, K_A_O, K_A_HCO, b_A_O, b_A_NO } = this.kin_params;
-    
+    const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params;
+
     // Hydrolysis
-    rates[0] = X_H == 0 ? 0 : k_H * this._monod(X_S / X_H, K_X) * X_H;
+    rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * this._monod(X_S / X_H, K_X) * X_H;
 
     // Heterotrophs
-    rates[1] = k_STO * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
-    rates[2] = k_STO * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_S, K_S) * X_H;
-    rates[3] = X_H == 0 ? 0 : mu_H_max * this._monod(S_O, K_O) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
-    rates[4] = X_H == 0 ? 0 : mu_H_max * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
-    rates[5] = b_H_O * this._monod(S_O, K_O) * X_H;
-    rates[6] = b_H_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
-    rates[7] = b_STO_O * this._monod(S_O, K_O) * X_H;
-    rates[8] = b_STO_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
+    rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
+    rates[2] = this._arrhenius(k_STO, theta_STO, T) * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_S, K_S) * X_H;
+    rates[3] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * this._monod(S_O, K_O) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
+    rates[4] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
+    rates[5] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
+    rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
+    rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
+    rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
     
     // Autotrophs
-    rates[9] = mu_A_max * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
-    rates[10] = b_A_O * this._monod(S_O, K_O) * X_A;
-    rates[11] = b_A_NO * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
+    rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
+    rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
+    rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
     
     return rates;
   }
 
   /**
-   * Computes the change in concentrations of reaction species based on the current state.
+   * Computes the change in concentrations of reaction species based on the current state and temperature.
    * @param {Array} state - State vector containing concentrations of reaction species.
+   * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
    * @returns {Array} - Change in reaction species concentrations.
    */
-  compute_dC(state) { // compute changes in concentrations
+  compute_dC(state, T = 20) { // compute changes in concentrations
     // state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
-    return math.multiply(this.stoi_matrix, this.compute_rates(state));
+    return math.multiply(this.stoi_matrix, this.compute_rates(state, T));
   }
 }