diff --git a/R/prj_chemical_system.R b/R/prj_chemical_system.R
index 264c865fcf486bcb2fc2ce7bbe85624e74f730e9..d563d33b2ca92b844e7cbf3e5e5b37f97033558e 100644
--- a/R/prj_chemical_system.R
+++ b/R/prj_chemical_system.R
@@ -14,6 +14,8 @@
 #' @param equilibrium_reactants Optional: list, prj_phase_component:
 #' @param surface Optional:
 #' @param user_punch Optional:
+#' @param linear_solver Optional:
+#' @param exchangers Optional: prj_exchangers
 #' @example man/examples/ex_prj_chemical_system.R
 #' @export
 prj_chemical_system <- function(chemical_solver,
@@ -25,7 +27,9 @@ prj_chemical_system <- function(chemical_solver,
                                    rates = NULL,
                                    equilibrium_reactants = NULL,
                                    surface = NULL,
-                                   user_punch = NULL) {
+                                   user_punch = NULL,
+                                   linear_solver = NULL,
+                                   exchangers = NULL) {
 
     # Add coercing utility here
 
@@ -38,7 +42,9 @@ prj_chemical_system <- function(chemical_solver,
                                rates,
                                equilibrium_reactants,
                                surface,
-                               user_punch)
+                               user_punch,
+                               linear_solver,
+                               exchangers)
 }
 
 
@@ -51,13 +57,17 @@ new_prj_chemical_system <- function(chemical_solver,
                                        rates = NULL,
                                        equilibrium_reactants = NULL,
                                        surface = NULL,
-                                       user_punch = NULL) {
+                                       user_punch = NULL,
+                                       linear_solver = NULL,
+                                       exchangers = NULL) {
 
 
     are_strings(chemical_solver,
                 database)
 
     assertthat::assert_that(class(solution) == "prj_solution")
+    assertthat::assert_that(is.null(exchangers) |
+                                class(exchangers) == "prj_exchangers")
 
     are_null_or_strings(mesh)
 
@@ -95,6 +105,8 @@ new_prj_chemical_system <- function(chemical_solver,
                    equilibrium_reactants = equilibrium_reactants,
                    surface = surface,
                    user_punch = user_punch,
+                   linear_solver = linear_solver,
+                   exchangers = exchangers,
                    xpath = "chemical_system",
                    attr_names = c("chemical_solver"),
                    flatten_on_exp = character()
@@ -141,7 +153,6 @@ new_prj_solution <- function(temperature,
                                 pe,
                                 components,
                                 charge_balance = NULL) {
-
     are_numbers(temperature,
                 pressure,
                 pe)
@@ -301,3 +312,27 @@ new_prj_rate <- function(kinetic_reactant,
         class = "prj_rate"
     )
 }
+
+#'prj_exchangers
+#'@description tag: exchangers
+#'@param exchange_site list
+#'@export
+prj_exchangers <- function(exchange_site) {
+
+    # Add coercing utility here
+
+    new_prj_exchangers(exchange_site)
+}
+
+new_prj_exchangers <- function(exchange_site) {
+
+    assertthat::assert_that(is.list(exchange_site))
+
+    structure(list(exchange_site = exchange_site,
+                   xpath = "chemical_system/exchangers",
+                   attr_names = character(),
+                   flatten_on_exp = character()
+    ),
+    class = "prj_exchangers"
+    )
+}
diff --git a/R/prj_medium.R b/R/prj_medium.R
index eee61a4225d51c7f3130d7ff2ea925e5541ff23b..ba966311b7d3f540f5598b505ee9e09c50c4ff39 100644
--- a/R/prj_medium.R
+++ b/R/prj_medium.R
@@ -79,8 +79,13 @@ new_prj_medium <- function(phases = NULL,
 #' @param intrinsic_permeability Optional:
 #' @param initial_aperture Optional:
 #' @param mean_frac_distance Optional:
+#' @param mean_frac_distances Optional: numeric vector
 #' @param threshold_strain Optional:
+#' @param threshold_strains Optional: numeric vector
 #' @param fracture_normal Optional:
+#' @param fracture_normals Optional: numeric vector
+#' @param fracture_rotation_xy Optional:
+#' @param fracture_rotation_yz Optional:
 #' @param reference_permeability Optional:
 #' @param fitting_factor Optional:
 #' @param cohesion Optional:
@@ -117,8 +122,13 @@ prj_pr_property <- function(name,
                             intrinsic_permeability = NULL,
                             initial_aperture = NULL,
                             mean_frac_distance = NULL,
+                            mean_frac_distances = NULL,
                             threshold_strain = NULL,
+                            threshold_strains = NULL,
                             fracture_normal = NULL,
+                            fracture_normals = NULL,
+                            fracture_rotation_xy = NULL,
+                            fracture_rotation_yz = NULL,
                             reference_permeability = NULL,
                             fitting_factor = NULL,
                             cohesion = NULL,
@@ -174,8 +184,13 @@ prj_pr_property <- function(name,
                         intrinsic_permeability,
                         initial_aperture,
                         mean_frac_distance,
+                        mean_frac_distances,
                         threshold_strain,
+                        threshold_strains,
                         fracture_normal,
+                        fracture_normals,
+                        fracture_rotation_xy,
+                        fracture_rotation_yz,
                         reference_permeability,
                         fitting_factor,
                         cohesion,
@@ -213,8 +228,13 @@ new_prj_pr_property <- function(name,
                                 intrinsic_permeability = NULL,
                                 initial_aperture = NULL,
                                 mean_frac_distance = NULL,
+                                mean_frac_distances = NULL,
                                 threshold_strain = NULL,
+                                threshold_strains = NULL,
                                 fracture_normal = NULL,
+                                fracture_normals = NULL,
+                                fracture_rotation_xy = NULL,
+                                fracture_rotation_yz = NULL,
                                 reference_permeability = NULL,
                                 fitting_factor = NULL,
                                 cohesion = NULL,
@@ -277,8 +297,13 @@ new_prj_pr_property <- function(name,
                    intrinsic_permeability = intrinsic_permeability,
                    initial_aperture = initial_aperture,
                    mean_frac_distance = mean_frac_distance,
+                   mean_frac_distances = mean_frac_distances,
                    threshold_strain = threshold_strain,
+                   threshold_strains = threshold_strains,
                    fracture_normal = fracture_normal,
+                   fracture_normals = fracture_normals,
+                   fracture_rotation_xy = fracture_rotation_xy,
+                   fracture_rotation_yz = fracture_rotation_yz,
                    reference_permeability = reference_permeability,
                    fitting_factor = fitting_factor,
                    cohesion = cohesion,
@@ -414,7 +439,10 @@ prj_ph_property <- function(name,
                             ...) {
 
     #Coerce input
-    value <- coerce_string_to_numeric(value)
+    if (!is.list(value)) {
+        value <- coerce_string_to_numeric(value)
+    }
+
     reference_value <- coerce_string_to_numeric(reference_value)
     exponents <- coerce_string_to_numeric(exponents)
     offset <- coerce_string_to_numeric(offset)
@@ -457,7 +485,11 @@ new_prj_ph_property <- function(name,
     are_strings(name,
                 type)
 
-    are_null_or_numeric(value)
+    if (is.list(value)) {
+        are_null_or_strings(value[[1]])
+    } else {
+        are_null_or_numeric(value)
+    }
 
     are_null_or_numbers(
         reference_value,
@@ -560,34 +592,79 @@ new_prj_component <- function(name,
 #' @param type string: Property type
 #' @param value Optional: string | double: ...
 #' @param parameter_name Optional:
+#' @param reference_diffusion Optional: character
+#' @param activation_energy Optional: string | double
+#' @param reference_temperature Optional: numeric
 #' @example man/examples/ex_prj_com_property.R
 #' @export
 prj_com_property <- function(name,
-                                type,
-                                value = NULL,
-                                parameter_name = NULL) {
+                             type,
+                             value = NULL,
+                             parameter_name = NULL,
+                             reference_diffusion = NULL,
+                             activation_energy = NULL,
+                             reference_temperature = NULL,
+                             triple_temperature = NULL,
+                             triple_pressure = NULL,
+                             critical_temperature = NULL,
+                             critical_pressure = NULL,
+                             reference_pressure= NULL) {
+
+
 
     #Coerce input
     value <- coerce_string_to_numeric(value)
+    activation_energy <- coerce_string_to_numeric(activation_energy)
+    reference_temperature <- coerce_string_to_numeric(reference_temperature)
+    triple_temperature <- coerce_string_to_numeric(triple_temperature)
+    triple_pressure <- coerce_string_to_numeric(triple_pressure)
+    critical_temperature <- coerce_string_to_numeric(critical_temperature)
+    critical_pressure <- coerce_string_to_numeric(critical_pressure)
+    reference_pressure <- coerce_string_to_numeric(reference_pressure)
 
     new_prj_com_property(name,
-                            type,
-                            value,
-                            parameter_name)
+                         type,
+                         value,
+                         parameter_name,
+                         reference_diffusion,
+                         activation_energy,
+                         reference_temperature,
+                         triple_temperature,
+                         triple_pressure,
+                         critical_temperature,
+                         critical_pressure,
+                         reference_pressure
+                         )
 }
 
 
 new_prj_com_property <- function(name,
-                                    type,
-                                    value = NULL,
-                                    parameter_name = NULL) {
+                                 type,
+                                 value = NULL,
+                                 parameter_name = NULL,
+                                 reference_diffusion = NULL,
+                                 activation_energy = NULL,
+                                 reference_temperature = NULL,
+                                 triple_temperature = NULL,
+                                 triple_pressure = NULL,
+                                 critical_temperature = NULL,
+                                 critical_pressure = NULL,
+                                 reference_pressure= NULL) {
 
 
     assertthat::assert_that(assertthat::is.string(name))
     assertthat::assert_that(assertthat::is.string(type))
 
     are_null_or_numbers(value)
+    are_null_or_numeric(activation_energy)
+    are_null_or_numeric(reference_temperature)
+    are_null_or_numeric(triple_temperature)
+    are_null_or_numeric(triple_pressure)
+    are_null_or_numeric(critical_temperature)
+    are_null_or_numeric(critical_pressure)
+    are_null_or_numeric(reference_pressure)
     are_null_or_strings(parameter_name)
+    are_null_or_strings(reference_diffusion)
 
     structure(
         list(
@@ -595,6 +672,14 @@ new_prj_com_property <- function(name,
             type = type,
             value = value,
             parameter_name = parameter_name,
+            reference_diffusion = reference_diffusion,
+            activation_energy = activation_energy,
+            reference_temperature = reference_temperature,
+            triple_temperature = triple_temperature,
+            triple_pressure = triple_pressure,
+            critical_temperature = critical_temperature,
+            critical_pressure = critical_pressure,
+            reference_pressure = reference_pressure,
             xpath = paste0("media/medium/phases/phase/components/component/",
                            "properties/property"),
             attr_names = character(),
diff --git a/R/prj_parameter.R b/R/prj_parameter.R
index 8092db5249358ff53221cdee2130e4146e99b586..5c23130a4eee315c05d34827c09fd62761d991ef 100644
--- a/R/prj_parameter.R
+++ b/R/prj_parameter.R
@@ -15,6 +15,8 @@
 #' @param mesh Optional: string:
 #' @param time_series Optional: list:
 #' @param use_local_coordinate_system Optional: string, "true" | "false":
+#' @param range Optional: list
+#' @param seed Optional: number
 #' @param ... Optional: for index_values and expression tags (since there can be
 #' multiple)
 #' @example man/examples/ex_prj_parameter.R
@@ -30,6 +32,8 @@ prj_parameter <- function(name,
                              mesh = NULL,
                              time_series = NULL,
                              use_local_coordinate_system = NULL,
+                             range = NULL,
+                             seed = NULL,
                              ...) {
 
     #Coerce input
@@ -41,6 +45,9 @@ prj_parameter <- function(name,
     index_values <- ellipsis_list[names(ellipsis_list) == "index_values"]
     expression <- ellipsis_list[names(ellipsis_list) == "expression"]
 
+    range <- coerce_string_to_numeric(range)
+    seed <- coerce_string_to_numeric(seed)
+
     new_prj_parameter(name,
                          type,
                          value,
@@ -53,7 +60,9 @@ prj_parameter <- function(name,
                          mesh,
                          expression,
                          time_series,
-                         use_local_coordinate_system)
+                         use_local_coordinate_system,
+                         range,
+                         seed)
 }
 
 
@@ -69,7 +78,9 @@ new_prj_parameter <- function(name,
                                  mesh = NULL,
                                  expression = NULL,
                                  time_series = NULL,
-                                 use_local_coordinate_system = NULL) {
+                                 use_local_coordinate_system = NULL,
+                                 range = NULL,
+                                 seed = NULL) {
 
     assertthat::assert_that(assertthat::is.string(name))
     assertthat::assert_that(assertthat::is.string(type))
@@ -104,6 +115,8 @@ new_prj_parameter <- function(name,
                                                       "parameter_name"))
         }
     }
+    are_null_or_numeric(range)
+    are_null_or_numbers(seed)
 
     structure(list(name = name,
                    type = type,
@@ -118,6 +131,8 @@ new_prj_parameter <- function(name,
                    expression = expression,
                    time_series = time_series,
                    use_local_coordinate_system = use_local_coordinate_system,
+                   range = range,
+                   seed = seed,
                    xpath = "parameters/parameter",
                    attr_names = character(),
                    flatten_on_exp = c("values"),
diff --git a/R/prj_process.R b/R/prj_process.R
index 264cfde1a206973b56de3219603fa02a465e2e62..14ac7a3f6cfb77274e5a8ecbd3e885d0ac46ae35 100644
--- a/R/prj_process.R
+++ b/R/prj_process.R
@@ -70,6 +70,14 @@
 #' @param pf_irrv Optional:
 #' @param micro_porosity Optional:
 #' @param explicit_hm_coupling_in_unsaturated_zone Optional:
+#' @param simplified_elasticity Optional: character
+#' @param chemically_induced_porosity_change Optional: character
+#' @param use_server_communication Optional:
+#' @param phasefield_model Optional:
+#' @param irreversible_threshold Optional: numeric
+#' @param tabular_file Optional: character
+#' @param temperature_field Optional: character
+#' @param use_stokes_brinkman_form Optional: character
 #' @param ... Optional: fracture_properties, constitutive_relation
 #' @example man/examples/ex_prj_process.R
 #' @export
@@ -137,8 +145,16 @@ prj_process <- function(name,
                            density_solid = NULL,
                            latent_heat_evaporation = NULL,
                            pf_irrv = NULL,
-                        micro_porosity = NULL,
+                           micro_porosity = NULL,
                            explicit_hm_coupling_in_unsaturated_zone = NULL,
+                           simplified_elasticity = NULL,
+                           chemically_induced_porosity_change = NULL,
+                           use_server_communication = NULL,
+                           phasefield_model = NULL,
+                           irreversible_threshold = NULL,
+                           tabular_file = NULL,
+                           temperature_field = NULL,
+                           use_stokes_brinkman_form = NULL,
                            ...){
 
     #Coerce input
@@ -171,6 +187,7 @@ prj_process <- function(name,
     split_method <- coerce_string_to_numeric(split_method)
     reg_param <- coerce_string_to_numeric(reg_param)
     pf_irrv <- coerce_string_to_numeric(pf_irrv)
+    irreversible_threshold <- coerce_string_to_numeric(irreversible_threshold)
 
     ellipsis_list <- list(...)
 
@@ -248,79 +265,95 @@ prj_process <- function(name,
         latent_heat_evaporation,
         pf_irrv,
         micro_porosity,
-        explicit_hm_coupling_in_unsaturated_zone
+        explicit_hm_coupling_in_unsaturated_zone,
+        simplified_elasticity,
+        chemically_induced_porosity_change,
+        use_server_communication,
+        phasefield_model,
+        irreversible_threshold,
+        tabular_file,
+        temperature_field,
+        use_stokes_brinkman_form
     )
 }
 
 
 new_prj_process <- function(name,
-                               type,
-                               integration_order,
-                               process_variables,
-                               secondary_variables = NULL,
-                               specific_body_force = NULL,
-                               constitutive_relation = NULL,
-                               solid_density = NULL,
-                               dimension = NULL,
-                               coupling_scheme = NULL,
-                               darcy_gravity = NULL,
-                               reference_temperature = NULL,
-                               fracture_model = NULL,
-                               fracture_properties = NULL,
-                               jacobian_assembler = NULL,
-                               internal_length = NULL,
-                               mass_lumping = NULL,
-                               porosity = NULL,
-                               calculatesurfaceflux = NULL,
-                               intrinsic_permeability = NULL,
-                               specific_storage = NULL,
-                               fluid_viscosity = NULL,
-                               biot_coefficient = NULL,
-                               fluid_density = NULL,
-                               initial_effective_stress = NULL,
-                               initial_fracture_effective_stress = NULL,
-                               phasefield_parameters = NULL,
-                               deactivate_matrix_in_flow = NULL,
-                               borehole_heat_exchangers = NULL,
-                               temperature = NULL,
-                               reactive_system = NULL,
-                               fluid_specific_heat_source = NULL,
-                               fluid_specific_isobaric_heat_capacity = NULL,
-                               solid_hydraulic_permeability = NULL,
-                               solid_specific_heat_source = NULL,
-                               solid_heat_conductivity = NULL,
-                               solid_specific_isobaric_heat_capacity = NULL,
-                               tortuosity = NULL,
-                               diffusion_coefficient = NULL,
-                               solid_density_dry = NULL,
-                               solid_density_initial = NULL,
-                               characteristic_pressure = NULL,
-                               characteristic_temperature = NULL,
-                               characteristic_vapour_mass_fraction = NULL,
-                               output_element_matrices = NULL,
-                               material_property = NULL,
-                               diffusion_coeff_component_b = NULL,
-                               diffusion_coeff_component_a = NULL,
-                               hydro_crack_scheme = NULL,
-                               at_num = NULL,
-                               initial_stress = NULL,
-                               split_method = NULL,
-                               reg_param = NULL,
-                               thermal_parameters = NULL,
-                               non_advective_form = NULL,
-                               fluid = NULL,
-                               porous_medium = NULL,
-                               decay_rate = NULL,
-                               fluid_reference_density = NULL,
-                               retardation_factor = NULL,
-                               solute_dispersivity_longitudinal = NULL,
-                               solute_dispersivity_transverse = NULL,
-                               molecular_diffusion_coefficient = NULL,
-                               density_solid = NULL,
-                               latent_heat_evaporation = NULL,
-                               pf_irrv = NULL,
+                            type,
+                            integration_order,
+                            process_variables,
+                            secondary_variables = NULL,
+                            specific_body_force = NULL,
+                            constitutive_relation = NULL,
+                            solid_density = NULL,
+                            dimension = NULL,
+                            coupling_scheme = NULL,
+                            darcy_gravity = NULL,
+                            reference_temperature = NULL,
+                            fracture_model = NULL,
+                            fracture_properties = NULL,
+                            jacobian_assembler = NULL,
+                            internal_length = NULL,
+                            mass_lumping = NULL,
+                            porosity = NULL,
+                            calculatesurfaceflux = NULL,
+                            intrinsic_permeability = NULL,
+                            specific_storage = NULL,
+                            fluid_viscosity = NULL,
+                            biot_coefficient = NULL,
+                            fluid_density = NULL,
+                            initial_effective_stress = NULL,
+                            initial_fracture_effective_stress = NULL,
+                            phasefield_parameters = NULL,
+                            deactivate_matrix_in_flow = NULL,
+                            borehole_heat_exchangers = NULL,
+                            temperature = NULL,
+                            reactive_system = NULL,
+                            fluid_specific_heat_source = NULL,
+                            fluid_specific_isobaric_heat_capacity = NULL,
+                            solid_hydraulic_permeability = NULL,
+                            solid_specific_heat_source = NULL,
+                            solid_heat_conductivity = NULL,
+                            solid_specific_isobaric_heat_capacity = NULL,
+                            tortuosity = NULL,
+                            diffusion_coefficient = NULL,
+                            solid_density_dry = NULL,
+                            solid_density_initial = NULL,
+                            characteristic_pressure = NULL,
+                            characteristic_temperature = NULL,
+                            characteristic_vapour_mass_fraction = NULL,
+                            output_element_matrices = NULL,
+                            material_property = NULL,
+                            diffusion_coeff_component_b = NULL,
+                            diffusion_coeff_component_a = NULL,
+                            hydro_crack_scheme = NULL,
+                            at_num = NULL,
+                            initial_stress = NULL,
+                            split_method = NULL,
+                            reg_param = NULL,
+                            thermal_parameters = NULL,
+                            non_advective_form = NULL,
+                            fluid = NULL,
+                            porous_medium = NULL,
+                            decay_rate = NULL,
+                            fluid_reference_density = NULL,
+                            retardation_factor = NULL,
+                            solute_dispersivity_longitudinal = NULL,
+                            solute_dispersivity_transverse = NULL,
+                            molecular_diffusion_coefficient = NULL,
+                            density_solid = NULL,
+                            latent_heat_evaporation = NULL,
+                            pf_irrv = NULL,
                             micro_porosity = NULL,
-                            explicit_hm_coupling_in_unsaturated_zone = NULL) {
+                            explicit_hm_coupling_in_unsaturated_zone = NULL,
+                            simplified_elasticity = NULL,
+                            chemically_induced_porosity_change = NULL,
+                            use_server_communication = NULL,
+                            phasefield_model = NULL,
+                            irreversible_threshold = NULL,
+                            tabular_file = NULL,
+                            temperature_field = NULL,
+                            use_stokes_brinkman_form = NULL) {
 
     #Basic validation
     assertthat::assert_that(assertthat::is.string(name))
@@ -416,7 +449,13 @@ new_prj_process <- function(name,
                                porosity,
                                deactivate_matrix_in_flow,
                                output_element_matrices,
-                               reference_temperature)
+                               reference_temperature,
+                               simplified_elasticity,
+                               chemically_induced_porosity_change,
+                               use_server_communication,
+                               tabular_file,
+                               temperature_field,
+                               use_stokes_brinkman_form)
 
 
     are_null_or_numbers(dimension,
@@ -437,7 +476,8 @@ new_prj_process <- function(name,
                                at_num,
                                split_method,
                                reg_param,
-                               pf_irrv)
+                               pf_irrv,
+                               irreversible_threshold)
 
     if(!is.null(mass_lumping)){
         mass_lumping <- stringr::str_remove_all(mass_lumping, "[:space:]*")
@@ -521,6 +561,15 @@ new_prj_process <- function(name,
             micro_porosity = micro_porosity,
             explicit_hm_coupling_in_unsaturated_zone =
                 explicit_hm_coupling_in_unsaturated_zone,
+            simplified_elasticity = simplified_elasticity,
+            chemically_induced_porosity_change =
+                chemically_induced_porosity_change,
+            use_server_communication = use_server_communication,
+            phasefield_model = phasefield_model,
+            irreversible_threshold = irreversible_threshold,
+            tabular_file = tabular_file,
+            temperature_field = temperature_field,
+            use_stokes_brinkman_form = use_stokes_brinkman_form,
             xpath = "processes/process",
             attr_names = c("secondary_variable"),
             flatten_on_exp = c("specific_body_force"),
@@ -985,7 +1034,7 @@ new_prj_jacobian_assembler <- function(type,
 prj_phasefield_parameters <- function(residual_stiffness,
                                          crack_resistance,
                                          crack_length_scale,
-                                         kinetic_coefficient,
+                                         kinetic_coefficient = NULL,
                                          history_field = NULL) {
 
     # Add coercing utility here
@@ -1001,15 +1050,15 @@ prj_phasefield_parameters <- function(residual_stiffness,
 new_prj_phasefield_parameters <- function(residual_stiffness,
                                              crack_resistance,
                                              crack_length_scale,
-                                             kinetic_coefficient,
+                                             kinetic_coefficient = NULL,
                                              history_field = NULL) {
 
     are_strings(residual_stiffness,
                        crack_resistance,
-                       crack_length_scale,
-                       kinetic_coefficient)
+                       crack_length_scale)
 
-    are_null_or_strings(history_field)
+    are_null_or_strings(kinetic_coefficient,
+                        history_field)
 
     structure(list(residual_stiffness = residual_stiffness,
                    crack_resistance = crack_resistance,
@@ -1094,6 +1143,7 @@ validate_process_variables <- function(process_variables){
                         "gas_pressure",
                         "capillary_pressure",
                         "liquid_pressure",
+                        "liquid_velocity",
                         "overall_mass_density")
 
     for(i in seq_len(length(process_variables))){
diff --git a/R/prj_time_loop.R b/R/prj_time_loop.R
index 17a728e60c2255e14857a3ec944dcb10596e535a..5346cba279dd9b87b7874cfb3f56dfc6d94f1f92 100644
--- a/R/prj_time_loop.R
+++ b/R/prj_time_loop.R
@@ -141,6 +141,7 @@ new_prj_tl_process <- function(ref,
 #' @param output_iteration_results Optional: string: Either "true" or "false"
 #' @param meshes Optional: character: A vector of mesh names
 #' @param fixed_output_times Optional: string | numeric:
+#' @param hdf Optional: numeric
 #' @example man/examples/ex_prj_output.R
 #' @export
 prj_output <- function(type,
@@ -152,7 +153,8 @@ prj_output <- function(type,
                           data_mode = NULL,
                           output_iteration_results = NULL,
                           meshes = NULL,
-                          fixed_output_times = NULL) {
+                          fixed_output_times = NULL,
+                          hdf = NULL) {
 
     #Coerce input
     fixed_output_times <- coerce_string_to_numeric(fixed_output_times)
@@ -170,7 +172,8 @@ prj_output <- function(type,
                          data_mode,
                          output_iteration_results,
                          meshes,
-                         fixed_output_times)
+                         fixed_output_times,
+                         hdf)
 }
 
 
@@ -183,7 +186,8 @@ new_prj_output <- function(type,
                                  data_mode = NULL,
                                  output_iteration_results = NULL,
                                  meshes = NULL,
-                                 fixed_output_times = NULL) {
+                                 fixed_output_times = NULL,
+                                 hdf = NULL) {
 
     assertthat::assert_that(assertthat::is.string(type))
     assertthat::assert_that(assertthat::is.string(prefix))
@@ -218,6 +222,7 @@ new_prj_output <- function(type,
              output_iteration_results = output_iteration_results,
              meshes = meshes,
              fixed_output_times = fixed_output_times,
+             hdf = hdf,
              xpath = "time_loop/output",
              attr_names = character(),
              flatten_on_exp = c("fixed_output_times")
diff --git a/R/sim_utils.R b/R/sim_utils.R
index 3b0e424a05882ae7db92bf1bfe4fc5c38d40de83..3085e99379d0d48a1161dc60d09f8d91372baf28 100644
--- a/R/sim_utils.R
+++ b/R/sim_utils.R
@@ -374,7 +374,7 @@ run_all_benchmarks <- function(path,
 
 
     # Filter nonexisting files from prj_paths
-    nonexisting_prj_paths <- prj_paths[!file.exists(prj_paths)]
+    nonexisting_prj_paths <- prj_paths[!file.exists(unlist(prj_paths))]
     prj_paths <- prj_paths[!prj_paths %in% nonexisting_prj_paths]
 
 
diff --git a/data/xpaths_for_classes.rda b/data/xpaths_for_classes.rda
index f0fc9f1b5cbfad31485df45ec0b2dd445880ec11..bc7d290fafb42ee7f1b21dbdde5f385f4fbf65c5 100644
Binary files a/data/xpaths_for_classes.rda and b/data/xpaths_for_classes.rda differ