diff --git a/doc/RTD/source/HydroSchemes/hopkins_sph.rst b/doc/RTD/source/HydroSchemes/hopkins_sph.rst
index bcc51e0ad96b18956f1c8e54f7bf2bf3b352c138..e4f1479230df96eabaa1fe16994960059858613b 100644
--- a/doc/RTD/source/HydroSchemes/hopkins_sph.rst
+++ b/doc/RTD/source/HydroSchemes/hopkins_sph.rst
@@ -28,3 +28,9 @@ scheme it includes a Monaghan AV scheme and a Balsara switch.
 .. code-block:: bash
    
    ./configure --with-hydro="pressure-energy"
+
+Both of the above schemes use a very simple, fixed artificial viscosity, only
+the ``SPH:viscosity_alpha`` parameter has any effect for this scheme. This will
+change the strength of the artificial viscosity throughout the simulation, and
+has a default of 0.8.
+
diff --git a/doc/RTD/source/HydroSchemes/minimal_sph.rst b/doc/RTD/source/HydroSchemes/minimal_sph.rst
index 1a16a23360aaba8b28920150af0d4f4b05c74c2f..bbcbe026b56381c007f58920f31115f9f9160d05 100644
--- a/doc/RTD/source/HydroSchemes/minimal_sph.rst
+++ b/doc/RTD/source/HydroSchemes/minimal_sph.rst
@@ -10,11 +10,17 @@ Minimal (Density-Energy) SPH
    :caption: Contents:
 
 This scheme is a textbook implementation of Density-Energy SPH, and can be used
-as a pedagogical example. It also implements a Monaghan AV scheme, like the
-GADGET-2 scheme. It uses very similar equations, but differs in implementation
-details; namely it tracks the internal energy \(u\) as the thermodynamic
-variable, rather than entropy \(A\). To use the minimal scheme, use
+as a pedagogical example. It also implements a Monaghan AV scheme with a
+Balsara switch, like the GADGET-2 scheme. It uses very similar equations, but
+differs in implementation details; namely it tracks the internal energy \(u\)
+as the thermodynamic variable, rather than entropy \(A\). To use the minimal
+scheme, use
 
 .. code-block:: bash
 
     ./configure --with-hydro="minimal"
+
+As it uses a very simple, fixed artificial viscosity, only the
+``SPH:viscosity_alpha`` parameter has any effect for this scheme. This will
+change the strength of the artificial viscosity throughout the simulation,
+and has a default of 0.8.
diff --git a/doc/RTD/source/HydroSchemes/traditional_sph.rst b/doc/RTD/source/HydroSchemes/traditional_sph.rst
index c69ea5f60644119b8590414ffe00a75246de49a6..455e8bebe516bd9be9f6df889f1ead2088ca94d2 100644
--- a/doc/RTD/source/HydroSchemes/traditional_sph.rst
+++ b/doc/RTD/source/HydroSchemes/traditional_sph.rst
@@ -15,3 +15,8 @@ a Monaghan artificial viscosity scheme and Balsara switch.
 To use this hydro scheme, you need no extra configuration options -- it is the
 default!
 
+As it uses a very simple, fixed artificial viscosity, only the
+``SPH:viscosity_alpha`` parameter has any effect for this scheme. This will
+change the strength of the artificial viscosity throughout the simulation,
+and has a default of 0.8.
+
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index ab422fd8f2f29a6b9fd200d4cfd6ea26ca120497..e1773154328fdbcb2ced4efae3e4787257dd022f 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -34,6 +34,10 @@ SPH:
   minimal_temperature:   0        # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
   H_mass_fraction:       0.755    # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy. Default value is derived from the physical constants.
   H_ionization_temperature: 1e4   # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
+  viscosity_alpha:       0.8      # (Optional) Override for the initial value of the artificial viscosity. In schemes that have a fixed AV, this remains as alpha throughout the run.
+  viscosity_alpha_max:   2.0      # (Optional) Maximal value for the artificial viscosity in schemes that allow alpha to vary
+  viscosity_alpha_min:   0.1      # (Optional) Minimal value for the artificial viscosity in schemes that allow alpha to vary
+  viscosity_length:      0.1      # (Optional) Decay length for the artificial viscosity in schemes that allow alpha to vary
 
 # Parameters for the self-gravity scheme
 Gravity:
diff --git a/src/const.h b/src/const.h
index 821e671ac4ce61334b669f60a4785a7daf12d51d..e417b8ca3827ef87396706c56df36bb9bd3aed75 100644
--- a/src/const.h
+++ b/src/const.h
@@ -21,18 +21,11 @@
 #define SWIFT_CONST_H
 
 /* SPH Viscosity constants. */
-/* Cosmology defaults: a=0.8, b=3.0. Planetary defaults: a=1.5, b=4.0
+/* Cosmology default beta=3.0. Planetary default beta=4.0
+ * Alpha can be set in the parameter file.
  * Beta is defined as in e.g. Price (2010) Eqn (103) */
-#define const_viscosity_alpha 0.8f
 #define const_viscosity_beta 3.0f
 
-#define const_viscosity_alpha_min \
-  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
-#define const_viscosity_alpha_max \
-  2.0f /* Values taken from (Price,2004), not used in legacy gadget mode */
-#define const_viscosity_length \
-  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
-
 /* SPH Thermal conductivity constants. */
 #define const_conductivity_alpha \
   1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index b495af1d5b3ba1b91522bdeeac27fe77b262bd22..8c29acef0079830951e5ca365485b43214335e07 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -458,12 +458,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
  * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties. 
  * @param dt_alpha The time-step used to evolve non-cosmological quantities such
  *                 as the artificial viscosity.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp,
-    const struct cosmology *cosmo, const float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
   const float fac_mu = cosmo->a_factor_mu;
 
   /* Some smoothing length multiples. */
@@ -492,6 +494,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->force.balsara =
       normDiv_v / (normDiv_v + normRot_v + 0.0001f * fac_mu * fc * h_inv);
 
+  /* Set the AV property */
+  p->alpha = hydro_props->viscosity.alpha;
+
   /* Viscosity parameter decay time */
   /* const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
    */
@@ -500,8 +505,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   /* const float S = max(-normDiv_v, 0.f); */
 
   /* Compute the particle's viscosity parameter time derivative */
-  /* const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau + */
-  /*                         (const_viscosity_alpha_max - p->alpha) * S; */
+  /* const float alpha_dot = (hydro_props->viscosity.alpha_max) - p->alpha) / tau + */
+  /*                         (hydro_props->viscosity.alpha_max - p->alpha) * S; */
 
   /* Update particle's viscosity paramter */
   /* p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase; */  // MATTHIEU
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index 51c448e7a6fa56216dc49197d8b9240ee5acf5f1..69919c202223fdecc197a87178e59767c02ee16e 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -183,13 +183,6 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
       h_grpsph, "Viscosity Model",
       "Morris & Monaghan (1997), Rosswog, Davies, Thielemann & "
       "Piran (2000) with additional Balsara (1995) switch");
-  io_write_attribute_f(h_grpsph, "Viscosity alpha_min",
-                       const_viscosity_alpha_min);
-  io_write_attribute_f(h_grpsph, "Viscosity alpha_max",
-                       const_viscosity_alpha_max);
-  io_write_attribute_f(h_grpsph, "Viscosity beta", 2.f);
-  io_write_attribute_f(h_grpsph, "Viscosity decay length",
-                       const_viscosity_length);
 
   /* Time integration properties */
   io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 36db326cf9624a77e3f2b4d06b5f67542f1edc9e..16bfeb053dc030c0ac34622deab5b911a2fec197 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -472,12 +472,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
  * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties. 
  * @param dt_alpha The time-step used to evolve non-cosmological quantities such
  *                 as the artificial viscosity.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp,
-    const struct cosmology *cosmo, const float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
   const float fac_mu = cosmo->a_factor_mu;
 
@@ -502,8 +504,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float P_over_rho2 = pressure * rho_inv * rho_inv;
 
   /* Compute the Balsara switch */
+  /* Pre-multiply in the AV factor; hydro_props are not passed to the iact functions */
   const float balsara =
-      abs_div_v / (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
 
   /* Compute the "grad h" term */
   const float omega_inv =
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 05b74b6a88ca49969f14965dd55fee827627c4db..ccab5cffadc40b2467953a319aa7e3c73ef577c0 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -497,8 +497,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -620,8 +619,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -654,8 +652,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 }
 
 #ifdef WITH_VECTORIZATION
-static const vector const_viscosity_alpha_fac =
-    FILL_VEC(-0.25f * const_viscosity_alpha);
 
 /**
  * @brief Force interaction computed using 1 vector
@@ -746,7 +742,7 @@ runner_iact_nonsym_1_vec_force(
 
   /* Now construct the full viscosity term */
   rho_ij.v = vec_mul(vec_set1(0.5f), vec_add(pirho.v, pjrho.v));
-  visc.v = vec_div(vec_mul(const_viscosity_alpha_fac.v,
+  visc.v = vec_div(vec_mul(vec_set1(-0.25f),
                            vec_mul(v_sig.v, vec_mul(mu_ij.v, balsara.v))),
                    rho_ij.v);
 
@@ -937,11 +933,11 @@ runner_iact_nonsym_2_vec_force(
   rho_ij.v = vec_mul(vec_set1(0.5f), vec_add(pirho.v, pjrho.v));
   rho_ij_2.v = vec_mul(vec_set1(0.5f), vec_add(pirho.v, pjrho_2.v));
 
-  visc.v = vec_div(vec_mul(const_viscosity_alpha_fac.v,
+  visc.v = vec_div(vec_mul(vec_set1(-0.25f),
                            vec_mul(v_sig.v, vec_mul(mu_ij.v, balsara.v))),
                    rho_ij.v);
   visc_2.v =
-      vec_div(vec_mul(const_viscosity_alpha_fac.v,
+      vec_div(vec_mul(vec_set1(-0.25f),
                       vec_mul(v_sig_2.v, vec_mul(mu_ij_2.v, balsara_2.v))),
               rho_ij_2.v);
 
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 27aab4b897395615dcc7d0b3701e90604bd9538f..ec7d34f7ad8697f1d639ea4951011ddb06ec8833 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -200,8 +200,6 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
   io_write_attribute_s(
       h_grpsph, "Viscosity Model",
       "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
-  io_write_attribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
-  io_write_attribute_f(h_grpsph, "Viscosity beta", 3.f);
 }
 
 /**
diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index 4268a8b1a208bb1d08dfe0a838d3b59eb8cd120c..bc7c7ab83358873de92e0eea49c01a03573ae3dd 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -451,12 +451,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
  * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties. 
  * @param dt_alpha The time-step used to evolve non-cosmological quantities such
  *                 as the artificial viscosity.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part* restrict p, struct xpart* restrict xp,
-    const struct cosmology* cosmo, const float dt_alpha) {
+    struct part *restrict p, struct xpart *restrict xp,
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
   /* Initialise values that are used in the force loop */
   p->flux.momentum[0] = 0.0f;
diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h
index 392fb10a614f1d45227f64a459bd8d574cf6ead2..327794a4d72d3f635866225af06783548fe3117e 100644
--- a/src/hydro/GizmoMFV/hydro.h
+++ b/src/hydro/GizmoMFV/hydro.h
@@ -476,12 +476,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
  * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties. 
  * @param dt_alpha The time-step used to evolve non-cosmological quantities such
  *                 as the artificial viscosity.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part* restrict p, struct xpart* restrict xp,
-    const struct cosmology* cosmo, const float dt_alpha) {
+    struct part *restrict p, struct xpart *restrict xp,
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
   /* Initialise values that are used in the force loop */
   p->gravity.mflux[0] = 0.0f;
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index cfaf482278fb3624649832b2f315a039f947c7e2..8e1d77fe251f87f95701aa4d890cd7d243d8fa52 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -389,6 +389,10 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
   p->density.rho_dh = 0.f;
+  p->density.div_v = 0.f;
+  p->density.rot_v[0] = 0.f;
+  p->density.rot_v[1] = 0.f;
+  p->density.rot_v[2] = 0.f;
 }
 
 /**
@@ -424,6 +428,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   p->density.rho_dh *= h_inv_dim_plus_one;
   p->density.wcount *= h_inv_dim;
   p->density.wcount_dh *= h_inv_dim_plus_one;
+
+  const float rho_inv = 1.f / p->rho;
+  const float a_inv2 = cosmo->a2_inv;
+
+  /* Finish calculation of the (physical) velocity curl components */
+  p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+  p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+  p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
+
+  /* Finish calculation of the (physical) velocity divergence */
+  p->density.div_v *= h_inv_dim_plus_one * a_inv2 * rho_inv;
 }
 
 /**
@@ -451,6 +466,10 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
   p->density.wcount = kernel_root * h_inv_dim;
   p->density.rho_dh = 0.f;
   p->density.wcount_dh = 0.f;
+  p->density.div_v = 0.f;
+  p->density.rot_v[0] = 0.f;
+  p->density.rot_v[1] = 0.f;
+  p->density.rot_v[2] = 0.f;
 }
 
 /**
@@ -466,12 +485,24 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
  * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties. 
  * @param dt_alpha The time-step used to evolve non-cosmological quantities such
  *                 as the artificial viscosity.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp,
-    const struct cosmology *cosmo, const float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
+
+  const float fac_mu = cosmo->a_factor_mu;
+
+  /* Compute the norm of the curl */
+  const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
+                             p->density.rot_v[1] * p->density.rot_v[1] +
+                             p->density.rot_v[2] * p->density.rot_v[2]);
+
+  /* Compute the norm of div v */
+  const float abs_div_v = fabsf(p->density.div_v);
 
   /* Compute the pressure */
   const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
@@ -484,10 +515,17 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float grad_h_term =
       1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
 
+  /* Compute the Balsara switch */
+  /* Pre-multiply in the AV factor; hydro_props are not passed to the iact functions */
+  const float balsara =
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
+
   /* Update variables. */
   p->force.f = grad_h_term;
   p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
+  p->force.balsara = balsara;
 }
 
 /**
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index e9ad55516b09d1b28988e68f7316edbc7352be71..13f3e5cce65c976ea4cdd5ccaeda20debf260387 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -80,6 +80,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
   pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
   pj->density.wcount += wj;
   pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+
+  /* Compute dv dot r */
+  float dv[3], curlvr[3];
+
+  const float faci = mj * wi_dx * r_inv;
+  const float facj = mi * wj_dx * r_inv;
+
+  dv[0] = pi->v[0] - pj->v[0];
+  dv[1] = pi->v[1] - pj->v[1];
+  dv[2] = pi->v[2] - pj->v[2];
+  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+  pi->density.div_v -= faci * dvdr;
+  pj->density.div_v -= facj * dvdr;
+
+  /* Compute dv cross r */
+  curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+  curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+  curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+  pi->density.rot_v[0] += faci * curlvr[0];
+  pi->density.rot_v[1] += faci * curlvr[1];
+  pi->density.rot_v[2] += faci * curlvr[2];
+
+  pj->density.rot_v[0] += facj * curlvr[0];
+  pj->density.rot_v[1] += facj * curlvr[1];
+  pj->density.rot_v[2] += facj * curlvr[2];
 }
 
 /**
@@ -115,6 +142,27 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
   pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
   pi->density.wcount += wi;
   pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+  /* Compute dv dot r */
+  float dv[3], curlvr[3];
+
+  const float faci = mj * wi_dx * r_inv;
+
+  dv[0] = pi->v[0] - pj->v[0];
+  dv[1] = pi->v[1] - pj->v[1];
+  dv[2] = pi->v[2] - pj->v[2];
+  const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+  pi->density.div_v -= faci * dvdr;
+
+  /* Compute dv cross r */
+  curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+  curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+  curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+  pi->density.rot_v[0] += faci * curlvr[0];
+  pi->density.rot_v[1] += faci * curlvr[1];
+  pi->density.rot_v[2] += faci * curlvr[2];
 }
 
 /**
@@ -186,9 +234,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float cj = pj->force.soundspeed;
   const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
 
+  /* Grab balsara switches */
+  const float balsara_i = pi->force.balsara;
+  const float balsara_j = pj->force.balsara;
+
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+  const float visc = -0.25f * v_sig * (balsara_i + balsara_j) * mu_ij / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -302,9 +354,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float cj = pj->force.soundspeed;
   const float v_sig = ci + cj - const_viscosity_beta * mu_ij;
 
+  /* Grab balsara switches */
+  const float balsara_i = pi->force.balsara;
+  const float balsara_j = pj->force.balsara;
+
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+  const float visc = -0.25f * v_sig * (balsara_i + balsara_j) * mu_ij / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index c33f1b9a214cf9839f1acb965b686d4a4962865c..1d14a94f2d91bf259df54c875a32bf3072ad33b6 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -124,6 +124,12 @@ struct part {
       /*! Derivative of density with respect to h */
       float rho_dh;
 
+      /*! Velocity divergence */
+      float div_v;
+
+      /*! Velocity curl */
+      float rot_v[3];
+
     } density;
 
     /**
@@ -150,6 +156,9 @@ struct part {
       /*! Time derivative of smoothing length  */
       float h_dt;
 
+      /*! Balsara switch */
+      float balsara;
+
     } force;
   };
 
diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h
index 1001dbda8b6b80a9e1c621a724350be4ec9a55c9..84b91057ae905c120c1210dba6e180344d852dfa 100644
--- a/src/hydro/Planetary/hydro.h
+++ b/src/hydro/Planetary/hydro.h
@@ -488,14 +488,15 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
  * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties. 
  * @param dt_alpha The time-step used to evolve non-cosmological quantities such
  *                 as the artificial viscosity.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp,
-    const struct cosmology *cosmo, const float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
-#ifndef PLANETARY_SPH_NO_BALSARA
   const float fac_mu = cosmo->a_factor_mu;
 
   /* Compute the norm of the curl */
@@ -505,7 +506,6 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Compute the norm of div v */
   const float abs_div_v = fabsf(p->density.div_v);
-#endif
 
   /* Compute the pressure */
   const float pressure =
@@ -527,20 +527,20 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     grad_h_term = 0.f;
   }
 
-#ifndef PLANETARY_SPH_NO_BALSARA
   /* Compute the Balsara switch */
+#ifdef PLANETARY_SPH_NO_BALSARA
+  const float balsara = hydro_props->viscosity.alpha;
+#else
   const float balsara =
-      abs_div_v / (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h);
 #endif
 
   /* Update variables. */
   p->force.f = grad_h_term;
   p->force.pressure = pressure;
   p->force.soundspeed = soundspeed;
-
-#ifndef PLANETARY_SPH_NO_BALSARA
   p->force.balsara = balsara;
-#endif
 }
 
 /**
diff --git a/src/hydro/Planetary/hydro_iact.h b/src/hydro/Planetary/hydro_iact.h
index 149ad9dda2c2b63cc1cdca1e55280c4c7676b8b6..c76512265fbb240d74c7570fc8b6783d72e05db1 100644
--- a/src/hydro/Planetary/hydro_iact.h
+++ b/src/hydro/Planetary/hydro_iact.h
@@ -176,11 +176,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
                      (pi->v[1] - pj->v[1]) * dx[1] +
                      (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
 
-#ifndef PLANETARY_SPH_NO_BALSARA
   /* Balsara term */
   const float balsara_i = pi->force.balsara;
   const float balsara_j = pj->force.balsara;
-#endif
 
   /* Are the particles moving towards each other? */
   const float omega_ij = min(dvdr, 0.f);
@@ -193,12 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-#ifdef PLANETARY_SPH_NO_BALSARA
-  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
-#else
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
-#endif
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -300,11 +293,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
                      (pi->v[1] - pj->v[1]) * dx[1] +
                      (pi->v[2] - pj->v[2]) * dx[2] + a2_Hubble * r2;
 
-#ifndef PLANETARY_SPH_NO_BALSARA
   /* Balsara term */
   const float balsara_i = pi->force.balsara;
   const float balsara_j = pj->force.balsara;
-#endif
 
   /* Are the particles moving towards each other? */
   const float omega_ij = min(dvdr, 0.f);
@@ -319,12 +310,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-#ifdef PLANETARY_SPH_NO_BALSARA
-  const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
-#else
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
+  const float visc = -0.25f * v_sig * mu_ij *
                      (balsara_i + balsara_j) / rho_ij;
-#endif
 
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
diff --git a/src/hydro/Planetary/hydro_part.h b/src/hydro/Planetary/hydro_part.h
index a40d86a53be371f3d3f50cdb9d263732cb2bdf08..4087cef62e873231a556f82869a7f6d848c8d72c 100644
--- a/src/hydro/Planetary/hydro_part.h
+++ b/src/hydro/Planetary/hydro_part.h
@@ -126,13 +126,11 @@ struct part {
       /*! Derivative of density with respect to h */
       float rho_dh;
 
-#ifndef PLANETARY_SPH_NO_BALSARA
       /*! Velocity divergence. */
       float div_v;
 
       /*! Velocity curl. */
       float rot_v[3];
-#endif
 
     } density;
 
@@ -160,10 +158,8 @@ struct part {
       /*! Time derivative of smoothing length  */
       float h_dt;
 
-#ifndef PLANETARY_SPH_NO_BALSARA
       /*! Balsara switch */
       float balsara;
-#endif
 
     } force;
   };
diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h
index df7de2124e377ecfe41247429eba31a23575dc35..5cde4685a3194e9c5c6a88c1b43fc1cea1a98706 100644
--- a/src/hydro/PressureEnergy/hydro.h
+++ b/src/hydro/PressureEnergy/hydro.h
@@ -524,12 +524,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
  * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties. 
  * @param dt_alpha The time-step used to evolve non-cosmological quantities such
  *                 as the artificial viscosity.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp,
-    const struct cosmology *cosmo, const float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
   const float fac_mu = cosmo->a_factor_mu;
 
@@ -546,7 +548,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Compute the Balsara switch */
   const float balsara =
-      abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
 
   /* Compute the "grad h" term */
   const float common_factor = p->h / (hydro_dimension * p->density.wcount);
diff --git a/src/hydro/PressureEnergy/hydro_iact.h b/src/hydro/PressureEnergy/hydro_iact.h
index 78aa4e724d5b61a125a30bd900211ea284f0db94..2ed7fe8cb8112c42de933a2ca315966e67108f0a 100644
--- a/src/hydro/PressureEnergy/hydro_iact.h
+++ b/src/hydro/PressureEnergy/hydro_iact.h
@@ -252,8 +252,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
@@ -380,8 +379,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Convolve with the kernel */
   const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 051b93e81f012106eed173679d3723dd5058827d..980b264468601ece3e4a036b50243099a7282015 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -478,12 +478,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
  * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties. 
  * @param dt_alpha The time-step used to evolve non-cosmological quantities such
  *                 as the artificial viscosity.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp,
-    const struct cosmology *cosmo, const float dt_alpha) {
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
   const float fac_mu = cosmo->a_factor_mu;
 
@@ -503,7 +505,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Compute the Balsara switch */
   const float balsara =
-      abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
+      hydro_props->viscosity.alpha * abs_div_v /
+      (abs_div_v + curl_v + 0.0001f * soundspeed * fac_mu / p->h);
 
   /* Divide the pressure by the density squared to get the SPH term */
   const float rho_bar_inv = 1.f / p->rho_bar;
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index f93e99137b8ffc47c495669047d25ab17062c593..a018b39a99be5ed691485d93bd8dfd1735378bda 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -263,8 +263,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
@@ -377,8 +376,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   /* Now construct the full viscosity term */
   const float rho_ij = 0.5f * (rhoi + rhoj);
-  const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
-                     (balsara_i + balsara_j) / rho_ij;
+  const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
 
   /* Now, convolve with the kernel */
   const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index 3c54d0e6385480c5905aa769d51546827390951d..e9397bf6108b8bc16658157e424055274f05f23c 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -194,8 +194,6 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
   io_write_attribute_s(
       h_grpsph, "Viscosity Model",
       "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
-  io_write_attribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
-  io_write_attribute_f(h_grpsph, "Viscosity beta", 3.f);
 
   /* Time integration properties */
   io_write_attribute_f(h_grpsph, "Maximal Delta u change over dt",
diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h
index 0715223850023be90294e5ccaa43c197c8eab2d4..70d94278bc40e2382df673ecdb12117d931ff3dc 100644
--- a/src/hydro/Shadowswift/hydro.h
+++ b/src/hydro/Shadowswift/hydro.h
@@ -295,12 +295,14 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
  * @param p The particle to act upon
  * @param xp The extended particle data to act upon
  * @param cosmo The current cosmological model.
+ * @param hydro_props Hydrodynamic properties. 
  * @param dt_alpha The time-step used to evolve non-cosmological quantities such
  *                 as the artificial viscosity.
  */
 __attribute__((always_inline)) INLINE static void hydro_prepare_force(
-    struct part* restrict p, struct xpart* restrict xp,
-    const struct cosmology* cosmo, const float dt_alpha) {
+    struct part *restrict p, struct xpart *restrict xp,
+    const struct cosmology *cosmo, const struct hydro_props *hydro_props,
+    const float dt_alpha) {
 
   /* Initialize time step criterion variables */
   p->timestepvars.vmax = 0.0f;
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 534319b5effb6e10ce831c3504dfb0c0eb99d41f..5cd7ba60e4e7c2dff123366da5787f47cbc8d10b 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -40,6 +40,13 @@
 #define hydro_props_default_init_temp 0.f
 #define hydro_props_default_min_temp 0.f
 #define hydro_props_default_H_ionization_temperature 1e4
+#define hydro_props_default_viscosity_alpha 0.8f
+#define hydro_props_default_viscosity_alpha_min \
+  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
+#define hydro_props_default_viscosity_alpha_max \
+  2.0f /* Values taken from (Price,2004), not used in legacy gadget mode */
+#define hydro_props_default_viscosity_length \
+  0.1f /* Values taken from (Price,2004), not used in legacy gadget mode */
 
 /**
  * @brief Initialize the global properties of the hydro scheme.
@@ -112,6 +119,23 @@ void hydro_props_init(struct hydro_props *p,
   p->hydrogen_mass_fraction = parser_get_opt_param_double(
       params, "SPH:H_mass_fraction", default_H_fraction);
 
+  /* Read the artificial viscosity parameters from the file, if they exist */
+  p->viscosity.alpha = parser_get_opt_param_float(
+    params, "SPH:viscosity_alpha", hydro_props_default_viscosity_alpha
+  );
+  
+  p->viscosity.alpha_max = parser_get_opt_param_float(
+    params, "SPH:viscosity_alpha_max", hydro_props_default_viscosity_alpha_max
+  );
+
+  p->viscosity.alpha_min = parser_get_opt_param_float(
+    params, "SPH:viscosity_alpha_min", hydro_props_default_viscosity_alpha_min
+  );
+
+  p->viscosity.length = parser_get_opt_param_float(
+    params, "SPH:viscosity_length", hydro_props_default_viscosity_length
+  );
+
   /* Compute the initial energy (Note the temp. read is in internal units) */
   /* u_init = k_B T_init / (mu m_p (gamma - 1)) */
   double u_init = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
@@ -164,6 +188,11 @@ void hydro_props_print(const struct hydro_props *p) {
 
   message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
 
+  message("Artificial viscosity parameters set to alpha: %.3f, max: %.3f, "
+          "min: %.3f, length: %.3f.",
+          p->viscosity.alpha, p->viscosity.alpha_max,
+          p->viscosity.alpha_min, p->viscosity.length);
+
   message(
       "Hydrodynamic integration: Max change of volume: %.2f "
       "(max|dlog(h)/dt|=%f).",
@@ -182,6 +211,7 @@ void hydro_props_print(const struct hydro_props *p) {
   if (p->minimal_temperature != hydro_props_default_min_temp)
     message("Minimal gas temperature set to %f", p->minimal_temperature);
 
+
     // Matthieu: Temporary location for this i/o business.
 
 #ifdef PLANETARY_SPH
@@ -221,10 +251,53 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
                        p->hydrogen_mass_fraction);
   io_write_attribute_f(h_grpsph, "Hydrogen ionization transition temperature",
                        p->hydrogen_ionization_temperature);
-  io_write_attribute_f(h_grpsph, "Alpha viscosity", const_viscosity_alpha);
+  io_write_attribute_f(h_grpsph, "Alpha viscosity", p->viscosity.alpha);
+  io_write_attribute_f(h_grpsph, "Alpha viscosity (max)",
+                       p->viscosity.alpha_max);
+  io_write_attribute_f(h_grpsph, "Alpha viscosity (min)",
+                       p->viscosity.alpha_min);
+  io_write_attribute_f(h_grpsph, "Viscosity decay length", p->viscosity.length);
+  io_write_attribute_f(h_grpsph, "Beta viscosity", const_viscosity_beta);
 }
 #endif
 
+/**
+ * @brief Initialises a hydro_props struct with somewhat useful values for
+ *        the automated test suite. This is not intended for production use,
+ *        but rather to fill for the purposes of mocking.
+ * 
+ * @param p the struct
+ */
+void hydro_props_init_no_hydro(struct hydro_props *p) {
+  p->eta_neighbours = 1.2348;
+  p->h_tolerance = hydro_props_default_h_tolerance;
+  p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
+  const float delta_eta = p->eta_neighbours * (1.f + p->h_tolerance);
+  p->delta_neighbours =
+      (pow_dimension(delta_eta) - pow_dimension(p->eta_neighbours)) *
+      kernel_norm;
+  p->h_max = hydro_props_default_h_max;
+  p->max_smoothing_iterations = hydro_props_default_max_iterations;
+  p->CFL_condition = 0.1;
+  p->log_max_h_change = logf(powf(1.4, hydro_dimension_inv));
+  
+  /* These values are inconsistent and in a production run would probably lead
+     to a crash. Again, this function is intended for mocking use in unit tests
+     and is _not_ to be used otherwise! */
+  p->minimal_temperature = hydro_props_default_min_temp;
+  p->minimal_internal_energy = 0.f;
+  p->initial_temperature = hydro_props_default_init_temp;
+  p->initial_internal_energy = 0.f;
+
+  p->hydrogen_mass_fraction = 0.755;
+  p->hydrogen_ionization_temperature = hydro_props_default_H_ionization_temperature;
+
+  p->viscosity.alpha = hydro_props_default_viscosity_alpha;
+  p->viscosity.alpha_max= hydro_props_default_viscosity_alpha_max;
+  p->viscosity.alpha_min= hydro_props_default_viscosity_alpha_min;
+  p->viscosity.length= hydro_props_default_viscosity_length;
+}
+
 /**
  * @brief Write a hydro_props struct to the given FILE as a stream of bytes.
  *
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index cb3fba9724b6c2e01a616505345a99c2259a67a2..41686a67832e85c823d7aa30a30ffa6338108ee5 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -78,11 +78,27 @@ struct hydro_props {
   /*! Initial internal energy per unit mass */
   float initial_internal_energy;
 
-  /*! Primoridal hydrogen mass fraction for initial energy conversion */
+  /*! Primordial hydrogen mass fraction for initial energy conversion */
   float hydrogen_mass_fraction;
 
   /*! Temperature of the neutral to ionized transition of Hydrogen */
   float hydrogen_ionization_temperature;
+
+  /*! Artificial viscosity parameters */
+  struct {
+    /*! For the fixed, simple case. Also used to set the initial AV
+        coefficient for variable schemes. */
+    float alpha;
+
+    /*! Artificial viscosity (max) for the variable case (e.g. M&M) */
+    float alpha_max;
+
+    /*! Artificial viscosity (min) for the variable case (e.g. M&M) */
+    float alpha_min;
+
+    /*! The decay length of the artificial viscosity (used in M&M, etc.) */
+    float length;
+  } viscosity;
 };
 
 void hydro_props_print(const struct hydro_props *p);
@@ -99,4 +115,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p);
 void hydro_props_struct_dump(const struct hydro_props *p, FILE *stream);
 void hydro_props_struct_restore(const struct hydro_props *p, FILE *stream);
 
+/* Setup for tests */
+void hydro_props_init_no_hydro(struct hydro_props *p);
+
 #endif /* SWIFT_HYDRO_PROPERTIES */
diff --git a/src/runner.c b/src/runner.c
index 08f9e780b97e592ce7b53ae6bbbcede9b78f8926..33aaeb03ceebcae18440830882d75796e97234db 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -863,6 +863,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
   const int with_cosmology = (e->policy & engine_policy_cosmology);
   const double time_base = e->time_base;
   const struct cosmology *cosmo = e->cosmology;
+  const struct hydro_props *hydro_props = e->hydro_properties;
 
   TIMER_TIC;
 
@@ -901,7 +902,7 @@ void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer) {
         }
 
         /* Compute variables required for the force loop */
-        hydro_prepare_force(p, xp, cosmo, dt_alpha);
+        hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
 
         /* The particle force values are now set.  Do _NOT_
            try to read any particle density variables! */
@@ -1034,6 +1035,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             hydro_reset_gradient(p);
 
 #else
+            /* This needs to be extracted here because otherwise it is an
+             * undefined variable for the schemes that use the extra ghost. */
+            const struct hydro_props *hydro_props = e->hydro_properties;
+
             /* Calculate the time-step for passing to hydro_prepare_force, used
              * for the evolution of alpha factors (i.e. those involved in the
              * artificial viscosity and thermal conduction terms) */
@@ -1053,7 +1058,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
             /* As of here, particle force variables will be set. */
 
             /* Compute variables required for the force loop */
-            hydro_prepare_force(p, xp, cosmo, dt_alpha);
+            hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
 
             /* The particle force values are now set.  Do _NOT_
                try to read any particle density variables! */
@@ -1132,6 +1137,10 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
         hydro_reset_gradient(p);
 
 #else
+        /* This needs to be extracted here because otherwise it is an
+         * undefined variable for the schemes that use the extra ghost. */
+        const struct hydro_props *hydro_props = e->hydro_properties;
+
         /* Calculate the time-step for passing to hydro_prepare_force, used for
          * the evolution of alpha factors (i.e. those involved in the artificial
          * viscosity and thermal conduction terms) */
@@ -1150,7 +1159,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
         /* As of here, particle force variables will be set. */
 
         /* Compute variables required for the force loop */
-        hydro_prepare_force(p, xp, cosmo, dt_alpha);
+        hydro_prepare_force(p, xp, cosmo, hydro_props, dt_alpha);
 
         /* The particle force values are now set.  Do _NOT_
            try to read any particle density variables! */
diff --git a/tests/testActivePair.c b/tests/testActivePair.c
index a96e4f661c06b75c14fd3b96fc3d28376b0b959e..1990ea1e4531ac60ae8e8911f1142add6c66d265 100644
--- a/tests/testActivePair.c
+++ b/tests/testActivePair.c
@@ -33,7 +33,8 @@
 
 /* Typdef function pointer for interaction function. */
 typedef void (*interaction_func)(struct runner *, struct cell *, struct cell *);
-typedef void (*init_func)(struct cell *, const struct cosmology *);
+typedef void (*init_func)(struct cell *, const struct cosmology *,
+                          const struct hydro_props *);
 typedef void (*finalise_func)(struct cell *, const struct cosmology *);
 
 /**
@@ -170,7 +171,8 @@ void clean_up(struct cell *ci) {
  * @brief Initializes all particles field to be ready for a density calculation
  */
 void zero_particle_fields_density(struct cell *c,
-                                  const struct cosmology *cosmo) {
+                                  const struct cosmology *cosmo,
+                                  const struct hydro_props *hydro_props) {
   for (int pid = 0; pid < c->hydro.count; pid++) {
     hydro_init_part(&c->hydro.parts[pid], NULL);
   }
@@ -179,7 +181,8 @@ void zero_particle_fields_density(struct cell *c,
 /**
  * @brief Initializes all particles field to be ready for a force calculation
  */
-void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo) {
+void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo,
+                                const struct hydro_props *hydro_props) {
   for (int pid = 0; pid < c->hydro.count; pid++) {
     struct part *p = &c->hydro.parts[pid];
     struct xpart *xp = &c->hydro.xparts[pid];
@@ -219,7 +222,7 @@ void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo) {
 #endif /* PRESSURE-ENERGY */
 
     /* And prepare for a round of force tasks. */
-    hydro_prepare_force(p, xp, cosmo, 0.);
+    hydro_prepare_force(p, xp, cosmo, hydro_props, 0.);
     hydro_reset_acceleration(p);
   }
 }
@@ -299,8 +302,8 @@ void test_pair_interactions(struct runner *runner, struct cell **ci,
   runner_do_sort(runner, *cj, 0x1FFF, 0, 0);
 
   /* Zero the fields */
-  init(*ci, runner->e->cosmology);
-  init(*cj, runner->e->cosmology);
+  init(*ci, runner->e->cosmology, runner->e->hydro_properties);
+  init(*cj, runner->e->cosmology, runner->e->hydro_properties);
 
   /* Run the test */
   vec_interaction(runner, *ci, *cj);
@@ -315,8 +318,8 @@ void test_pair_interactions(struct runner *runner, struct cell **ci,
   /* Now perform a brute-force version for accuracy tests */
 
   /* Zero the fields */
-  init(*ci, runner->e->cosmology);
-  init(*cj, runner->e->cosmology);
+  init(*ci, runner->e->cosmology, runner->e->hydro_properties);
+  init(*cj, runner->e->cosmology, runner->e->hydro_properties);
 
   /* Run the brute-force test */
   serial_interaction(runner, *ci, *cj);
@@ -487,6 +490,7 @@ int main(int argc, char *argv[]) {
   struct space space;
   struct engine engine;
   struct cosmology cosmo;
+  struct hydro_props hydro_props;
   struct runner *runner;
   char c;
   static long long partId = 0;
@@ -571,6 +575,8 @@ int main(int argc, char *argv[]) {
 
   cosmology_init_no_cosmo(&cosmo);
   engine.cosmology = &cosmo;
+  hydro_props_init_no_hydro(&hydro_props);
+  engine.hydro_properties = &hydro_props;
 
   if (posix_memalign((void **)&runner, SWIFT_STRUCT_ALIGNMENT,
                      sizeof(struct runner)) != 0) {