diff --git a/src/Mod/Fem/femsolver/elmer/equations/flux.py b/src/Mod/Fem/femsolver/elmer/equations/flux.py
index 17edbbc575..4930a7c3cc 100644
--- a/src/Mod/Fem/femsolver/elmer/equations/flux.py
+++ b/src/Mod/Fem/femsolver/elmer/equations/flux.py
@@ -45,42 +45,30 @@ class Proxy(linear.Proxy, equationbase.FluxProxy):
 
     def __init__(self, obj):
         super(Proxy, self).__init__(obj)
-        obj.addProperty(
-            "App::PropertyBool",
-            "DiscontinuousGalerkin",
-            "Flux",
-            ""
-        )
+
         obj.addProperty(
             "App::PropertyBool",
             "AverageWithinMaterials",
             "Flux",
-            ""
+            "Enforces continuity within the same material\nin the 'Discontinuous Galerkin' discretization"
         )
         obj.addProperty(
             "App::PropertyBool",
             "CalculateFlux",
             "Flux",
-            ""
-        )
-        obj.CalculateFlux = True
-        obj.addProperty(
-            "App::PropertyString",
-            "FluxVariable",
-            "Flux",
-            "Insert variable name for flux calculation"
+            "Computes flux vector"
         )
         obj.addProperty(
             "App::PropertyBool",
             "CalculateFluxAbs",
             "Flux",
-            "Select calculation of abs of flux"
+            "Computes absolute of flux vector"
         )
         obj.addProperty(
             "App::PropertyBool",
             "CalculateFluxMagnitude",
             "Flux",
-            "Select calculation of magnitude of flux"
+            "Computes magnitude of flux vector field"
         )
         obj.addProperty(
             "App::PropertyBool",
@@ -92,21 +80,44 @@ class Proxy(linear.Proxy, equationbase.FluxProxy):
             "App::PropertyBool",
             "CalculateGradAbs",
             "Flux",
-            "Select calculation of abs of gradient"
+            "Computes absolute of gradient field"
         )
         obj.addProperty(
             "App::PropertyBool",
             "CalculateGradMagnitude",
             "Flux",
-            "Select calculation of magnitude of gradient"
+            "Computes magnitude of gradient field"
+        )
+        obj.addProperty(
+            "App::PropertyBool",
+            "DiscontinuousGalerkin",
+            "Flux",
+            "Enable if standard Galerkin approximation leads to\nunphysical results when there are discontinuities"
         )
         obj.addProperty(
             "App::PropertyBool",
             "EnforcePositiveMagnitude",
             "Flux",
-            "Select calculation of positive magnitude"
+            "If true, negative values of computed magnitude fields\nare a posteriori set to zero."
         )
+        obj.addProperty(
+            "App::PropertyString",
+            "FluxCoefficient",
+            "Flux",
+            "Name of proportionality coefficient\nto compute the flux"
+        )
+        obj.addProperty(
+            "App::PropertyString",
+            "FluxVariable",
+            "Flux",
+            "Variable name for flux calculation"
+        )
+        
         obj.Priority = 5
+        obj.CalculateFlux = True
+        # set defaults according to the Elmer manual
+        obj.FluxCoefficient = "Temperature"
+        obj.FluxVariable = "Heat Conductivity"
 
 
 class ViewProxy(linear.ViewProxy, equationbase.FluxViewProxy):
diff --git a/src/Mod/Fem/femsolver/elmer/writer.py b/src/Mod/Fem/femsolver/elmer/writer.py
index b649a09b7c..8257f882cb 100644
--- a/src/Mod/Fem/femsolver/elmer/writer.py
+++ b/src/Mod/Fem/femsolver/elmer/writer.py
@@ -523,16 +523,17 @@ class Writer(object):
                     activeIn = equation.References[0][1]
                 else:
                     activeIn = self._getAllBodies()
-                solverSection = self._getFlux(equation)
+                solverSection = self._getFluxSolver(equation)
                 for body in activeIn:
                     self._addSolver(body, solverSection)
 
-    def _getFlux(self, equation):
+    def _getFluxSolver(self, equation):
         s = self._createLinearSolver(equation)
+        # check if we need to update the equation
+        self._updateFluxSolver(equation)
+        # output the equation parameters
         s["Equation"] = "Flux Solver"  # equation.Name
         s["Procedure"] = sifio.FileAttr("FluxSolver/FluxSolver")
-        s["Flux Variable"] = equation.FluxVariable
-        s["Discontinuous Galerkin"] = equation.DiscontinuousGalerkin
         s["Average Within Materials"] = equation.AverageWithinMaterials
         s["Calculate Flux"] = equation.CalculateFlux
         s["Calculate Flux Abs"] = equation.CalculateFluxAbs
@@ -540,10 +541,75 @@ class Writer(object):
         s["Calculate Grad"] = equation.CalculateGrad
         s["Calculate Grad Abs"] = equation.CalculateGradAbs
         s["Calculate Grad Magnitude"] = equation.CalculateGradMagnitude
+        s["Discontinuous Galerkin"] = equation.DiscontinuousGalerkin
         s["Enforce Positive Magnitude"] = equation.EnforcePositiveMagnitude
+        s["Flux Coefficient"] = equation.FluxCoefficient
+        s["Flux Variable"] = equation.FluxVariable
         s["Stabilize"] = equation.Stabilize
         return s
 
+    def _updateFluxSolver(self, equation):
+        # updates older Flux equations
+        if not hasattr(equation, "AverageWithinMaterials"):
+            equation.addProperty(
+            "App::PropertyBool",
+            "AverageWithinMaterials",
+            "Flux",
+            "Enforces continuity within the same material\nin the 'Discontinuous Galerkin' discretization"
+            )
+        if hasattr(equation, "Bubbles"):
+            # Bubbles was removed because it is unused by Elmer for the flux solver
+            equation.removeProperty("Bubbles")
+        if not hasattr(equation, "CalculateFluxAbs"):
+            equation.addProperty(
+            "App::PropertyBool",
+            "CalculateFluxAbs",
+            "Flux",
+            "Computes absolute of flux vector"
+            )
+        if not hasattr(equation, "CalculateFluxMagnitude"):
+            equation.addProperty(
+            "App::PropertyBool",
+            "CalculateFluxMagnitude",
+            "Flux",
+            "Computes magnitude of flux vector field"
+            )
+        if not hasattr(equation, "CalculateGradAbs"):
+            equation.addProperty(
+            "App::PropertyBool",
+            "CalculateGradAbs",
+            "Flux",
+            "Computes absolute of gradient field"
+            )
+        if not hasattr(equation, "CalculateGradMagnitude"):
+            equation.addProperty(
+            "App::PropertyBool",
+            "CalculateGradMagnitude",
+            "Flux",
+            "Computes magnitude of gradient field"
+            )
+        if not hasattr(equation, "DiscontinuousGalerkin"):
+            equation.addProperty(
+            "App::PropertyBool",
+            "DiscontinuousGalerkin",
+            "Flux",
+            "Enable if standard Galerkin approximation leads to\nunphysical results when there are discontinuities"
+            )
+        if not hasattr(equation, "EnforcePositiveMagnitude"):
+            equation.addProperty(
+            "App::PropertyBool",
+            "EnforcePositiveMagnitude",
+            "Flux",
+            "If true, negative values of computed magnitude fields\nare a posteriori set to zero."
+            )
+        if not hasattr(equation, "FluxCoefficient"):
+            equation.addProperty(
+            "App::PropertyString",
+            "FluxCoefficient",
+            "Flux",
+            "Name of proportionality coefficient\nto compute the flux"
+            )
+
     def _handleElectricforce(self):
         activeIn = []
         for equation in self.solver.Group: