diff --git a/examples/fluid/falling_water_spheres_2d.jl b/examples/fluid/falling_water_spheres_2d.jl
index 68eae04574..667264d668 100644
--- a/examples/fluid/falling_water_spheres_2d.jl
+++ b/examples/fluid/falling_water_spheres_2d.jl
@@ -94,7 +94,7 @@ ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=1000)
 saving_callback = SolutionSavingCallback(dt=0.01, output_directory="out",
-                                         prefix="", write_meta_data=true)
+                                         prefix="")
 
 callbacks = CallbackSet(info_callback, saving_callback)
 
diff --git a/examples/fsi/falling_sphere_3d.jl b/examples/fsi/falling_sphere_3d.jl
index 53a709b243..0f49aa7b29 100644
--- a/examples/fsi/falling_sphere_3d.jl
+++ b/examples/fsi/falling_sphere_3d.jl
@@ -16,5 +16,4 @@ trixi_include(@__MODULE__,
               solid_smoothing_kernel=WendlandC2Kernel{3}(),
               sphere_type=RoundSphere(),
               output_directory="out", prefix="",
-              write_meta_data=false, # Files with meta data can't be read by meshio
               tspan=(0.0, 1.0), abstol=1e-6, reltol=1e-3)
diff --git a/examples/fsi/falling_spheres_2d.jl b/examples/fsi/falling_spheres_2d.jl
index 8c76b763eb..d59766befe 100644
--- a/examples/fsi/falling_spheres_2d.jl
+++ b/examples/fsi/falling_spheres_2d.jl
@@ -124,8 +124,7 @@ semi = Semidiscretization(fluid_system, boundary_system, solid_system_1, solid_s
 ode = semidiscretize(semi, tspan)
 
 info_callback = InfoCallback(interval=50)
-saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix="",
-                                         write_meta_data=true)
+saving_callback = SolutionSavingCallback(dt=0.02, output_directory="out", prefix="")
 
 callbacks = CallbackSet(info_callback, saving_callback)
 
diff --git a/examples/n_body/n_body_system.jl b/examples/n_body/n_body_system.jl
index 9493f2333f..b877059a41 100644
--- a/examples/n_body/n_body_system.jl
+++ b/examples/n_body/n_body_system.jl
@@ -105,7 +105,7 @@ end
 
 TrixiParticles.vtkname(system::NBodySystem) = "n-body"
 
-function TrixiParticles.write2vtk!(vtk, v, u, t, system::NBodySystem; write_meta_data=true)
+function TrixiParticles.write2vtk!(vtk, v, u, t, system::NBodySystem)
     (; mass) = system
 
     vtk["velocity"] = v
@@ -114,6 +114,10 @@ function TrixiParticles.write2vtk!(vtk, v, u, t, system::NBodySystem; write_meta
     return vtk
 end
 
+function TrixiParticles.add_system_data!(system_data, system::NBodySystem)
+    return system_data
+end
+
 function Base.show(io::IO, system::NBodySystem)
     print(io, "NBodySystem{", ndims(system), "}() with ")
     print(io, TrixiParticles.nparticles(system), " particles")
diff --git a/src/TrixiParticles.jl b/src/TrixiParticles.jl
index 87ce6d43f9..2dffe6a847 100644
--- a/src/TrixiParticles.jl
+++ b/src/TrixiParticles.jl
@@ -56,8 +56,8 @@ include("callbacks/callbacks.jl")
 # included separately. `gpu.jl` in turn depends on the semidiscretization type.
 include("general/semidiscretization.jl")
 include("general/gpu.jl")
-include("io/io.jl")
 include("preprocessing/preprocessing.jl")
+include("io/io.jl")
 include("visualization/recipes_plots.jl")
 
 export Semidiscretization, semidiscretize, restart_with!
diff --git a/src/callbacks/post_process.jl b/src/callbacks/post_process.jl
index 808b2eafc3..40b31998fa 100644
--- a/src/callbacks/post_process.jl
+++ b/src/callbacks/post_process.jl
@@ -266,7 +266,7 @@ function (pp::PostprocessCallback)(integrator)
 
         if isfinished(integrator) ||
            (pp.write_file_interval > 0 && backup_condition(pp, integrator))
-            write_postprocess_callback(pp)
+            write_postprocess_callback(pp, integrator)
         end
 
         # Tell OrdinaryDiffEq that `u` has not been modified
@@ -285,13 +285,14 @@ end
 end
 
 # After the simulation has finished, this function is called to write the data to a JSON file
-function write_postprocess_callback(pp::PostprocessCallback)
+function write_postprocess_callback(pp::PostprocessCallback, integrator)
     isempty(pp.data) && return
 
     mkpath(pp.output_directory)
 
     data = Dict{String, Any}()
-    write_meta_data!(data, pp.git_hash[])
+    data["meta"] = create_meta_data_dict(pp, integrator)
+
     prepare_series_data!(data, pp)
 
     time_stamp = ""
@@ -341,15 +342,6 @@ function create_series_dict(values, times, system_name="")
                 "time" => times)
 end
 
-function write_meta_data!(data, git_hash)
-    meta_data = Dict("solver_name" => "TrixiParticles.jl",
-                     "solver_version" => git_hash,
-                     "julia_version" => string(VERSION))
-
-    data["meta"] = meta_data
-    return data
-end
-
 function write_csv(abs_file_path, data)
     times = Float64[]
 
diff --git a/src/callbacks/solution_saving.jl b/src/callbacks/solution_saving.jl
index 39a0925210..0d017249c0 100644
--- a/src/callbacks/solution_saving.jl
+++ b/src/callbacks/solution_saving.jl
@@ -2,8 +2,7 @@
     SolutionSavingCallback(; interval::Integer=0, dt=0.0, save_times=Array{Float64, 1}([]),
                            save_initial_solution=true, save_final_solution=true,
                            output_directory="out", append_timestamp=false, prefix="",
-                           verbose=false, write_meta_data=true, max_coordinates=2^15,
-                           custom_quantities...)
+                           verbose=false, max_coordinates=2^15, custom_quantities...)
 
 
 Callback to save the current numerical solution in VTK format in regular intervals.
@@ -26,9 +25,8 @@ To ignore a custom quantity for a specific system, return `nothing`.
 - `save_final_solution=true`:   Save the final solution.
 - `output_directory="out"`:     Directory to save the VTK files.
 - `append_timestamp=false`:     Append current timestamp to the output directory.
-- 'prefix=""':                  Prefix added to the filename.
+- `prefix=""`:                  Prefix added to the filename.
 - `custom_quantities...`:       Additional user-defined quantities.
-- `write_meta_data=true`:       Write meta data.
 - `verbose=false`:              Print to standard IO when a file is written.
 - `max_coordinates=2^15`:       The coordinates of particles will be clipped if their
                                 absolute values exceed this threshold.
@@ -70,7 +68,6 @@ mutable struct SolutionSavingCallback{I, CQ}
     save_times            :: Vector{Float64}
     save_initial_solution :: Bool
     save_final_solution   :: Bool
-    write_meta_data       :: Bool
     verbose               :: Bool
     output_directory      :: String
     prefix                :: String
@@ -84,8 +81,9 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0,
                                 save_times=Float64[],
                                 save_initial_solution=true, save_final_solution=true,
                                 output_directory="out", append_timestamp=false,
-                                prefix="", verbose=false, write_meta_data=true,
-                                max_coordinates=Float64(2^15), custom_quantities...)
+                                prefix="", verbose=false,
+                                max_coordinates=Float64(2^15),
+                                custom_quantities...)
     if (dt > 0 && interval > 0) || (length(save_times) > 0 && (dt > 0 || interval > 0))
         throw(ArgumentError("Setting multiple save times for the same solution " *
                             "callback is not possible. Use either `dt`, `interval` or `save_times`."))
@@ -101,8 +99,8 @@ function SolutionSavingCallback(; interval::Integer=0, dt=0.0,
 
     solution_callback = SolutionSavingCallback(interval, Float64.(save_times),
                                                save_initial_solution, save_final_solution,
-                                               write_meta_data, verbose, output_directory,
-                                               prefix, max_coordinates, custom_quantities,
+                                               verbose, output_directory, prefix,
+                                               max_coordinates, custom_quantities,
                                                -1, Ref("UnknownVersion"))
 
     if length(save_times) > 0
@@ -133,6 +131,8 @@ function initialize_save_cb!(solution_callback::SolutionSavingCallback, u, t, in
     solution_callback.latest_saved_iter = -1
     solution_callback.git_hash[] = compute_git_hash()
 
+    write_meta_data(solution_callback, integrator)
+
     # Save initial solution
     if solution_callback.save_initial_solution
         solution_callback(integrator)
@@ -151,8 +151,8 @@ end
 
 # `affect!`
 function (solution_callback::SolutionSavingCallback)(integrator)
-    (; interval, output_directory, custom_quantities, write_meta_data, git_hash,
-     verbose, prefix, latest_saved_iter, max_coordinates) = solution_callback
+    (; interval, output_directory, custom_quantities, git_hash, verbose,
+     prefix, latest_saved_iter, max_coordinates) = solution_callback
 
     dvdu_ode = get_du(integrator)
     vu_ode = integrator.u
@@ -176,8 +176,8 @@ function (solution_callback::SolutionSavingCallback)(integrator)
 
     @trixi_timeit timer() "save solution" trixi2vtk(dvdu_ode, vu_ode, semi, integrator.t;
                                                     iter, output_directory, prefix,
-                                                    write_meta_data, git_hash=git_hash[],
-                                                    max_coordinates, custom_quantities...)
+                                                    git_hash=git_hash[], max_coordinates,
+                                                    custom_quantities...)
 
     # Tell OrdinaryDiffEq that `u` has not been modified
     u_modified!(integrator, false)
diff --git a/src/io/io.jl b/src/io/io.jl
index 380151bcff..680bbe9ac2 100644
--- a/src/io/io.jl
+++ b/src/io/io.jl
@@ -1,2 +1,321 @@
 include("write_vtk.jl")
 include("read_vtk.jl")
+
+# Handle "_" on optional prefix strings
+add_underscore_to_optional_prefix(str) = (str === "" ? "" : "$(str)_")
+# Same for optional postfix strings
+add_underscore_to_optional_postfix(str) = (str === "" ? "" : "_$(str)")
+
+function write_meta_data(callback::SolutionSavingCallback, integrator)
+    prefix = callback.prefix
+
+    meta_data = create_meta_data_dict(callback, integrator)
+
+    # write JSON file
+    output_directory = callback.output_directory
+    mkpath(output_directory)
+    json_file = joinpath(output_directory,
+                         "meta" * add_underscore_to_optional_postfix(prefix) * ".json")
+
+    open(json_file, "w") do file
+        JSON.print(file, meta_data, 2)
+    end
+end
+
+function create_meta_data_dict(callback, integrator)
+    git_hash = callback.git_hash
+    prefix = hasproperty(callback, :prefix) ? callback.prefix : ""
+    semi = integrator.p
+    names = system_names(semi.systems)
+
+    meta_data = Dict{String, Any}()
+
+    info = Dict{String, Any}()
+    add_simulation_info!(info, git_hash, integrator)
+    meta_data["simulation_info"] = info
+
+    systems = Dict{String, Any}()
+    foreach_system(semi) do system
+        idx = system_indices(system, semi)
+        name = add_underscore_to_optional_prefix(prefix) * names[idx]
+
+        system_data = Dict{String, Any}()
+        add_system_data!(system_data, system)
+
+        systems[name] = system_data
+    end
+    meta_data["system_data"] = systems
+
+    return meta_data
+end
+
+function add_simulation_info!(info, git_hash, integrator)
+    info["solver_name"] = "TrixiParticles.jl"
+    info["solver_version"] = git_hash[]
+    info["julia_version"] = string(VERSION)
+
+    info["time_integrator"] = Dict{String, Any}()
+    info["time_integrator"]["integrator_type"] = type2string(integrator.alg)
+    info["time_integrator"]["start_time"] = first(integrator.sol.prob.tspan)
+    info["time_integrator"]["final_time"] = last(integrator.sol.prob.tspan)
+    info["time_integrator"]["adaptive"] = integrator.opts.adaptive
+    if integrator.opts.adaptive
+        info["time_integrator"]["abstol"] = integrator.opts.abstol
+        info["time_integrator"]["reltol"] = integrator.opts.reltol
+        info["time_integrator"]["controller"] = type2string(integrator.opts.controller)
+    else
+        info["time_integrator"]["dt"] = integrator.dt
+        info["time_integrator"]["dt_max"] = integrator.opts.dtmax
+    end
+
+    info["technical_setup"] = Dict{String, Any}()
+    info["technical_setup"]["parallelization_backend"] = type2string(integrator.p.parallelization_backend)
+    info["technical_setup"]["number_of_threads"] = Threads.nthreads()
+end
+
+add_system_data!(system_data, data::Nothing) = system_data
+
+function add_system_data!(system_data, system::FluidSystem)
+    system_data["system_type"] = type2string(system)
+    system_data["particle_spacing"] = particle_spacing(system, 1)
+    system_data["density_calculator"] = type2string(system.density_calculator)
+    system_data["smoothing_kernel"] = type2string(system.smoothing_kernel)
+    system_data["smoothing_length"] = system.cache.smoothing_length
+    system_data["acceleration"] = system.acceleration
+    system_data["sound_speed"] = system_sound_speed(system)
+    system_data["pressure_acceleration_formulation"] = nameof(system.pressure_acceleration_formulation)
+    add_system_data!(system_data, shifting_technique(system))
+    add_system_data!(system_data, system.surface_tension)
+    add_system_data!(system_data, system.surface_normal_method)
+    add_system_data!(system_data, system.viscosity)
+    add_system_data!(system_data, system.correction)
+    add_system_data!(system_data, system_state_equation(system))
+    if hasfield(typeof(system), :density_diffusion)
+        add_system_data!(system_data, system.density_diffusion)
+    end
+    if hasfield(typeof(system), :alpha)
+        system_data["alpha"] = system.alpha
+    end
+end
+
+function add_system_data!(system_data, system::TotalLagrangianSPHSystem)
+    system_data["system_type"] = type2string(system)
+    system_data["particle_spacing"] = particle_spacing(system, 1)
+    system_data["smoothing_kernel"] = type2string(system.smoothing_kernel)
+    system_data["smoothing_length"] = system.smoothing_length
+    system_data["acceleration"] = system.acceleration
+    add_system_data!(system_data, system.boundary_model)
+    add_system_data!(system_data, system.viscosity)
+    add_system_data!(system_data, system.penalty_force)
+end
+
+function add_system_data!(system_data, system::BoundarySPHSystem)
+    system_data["system_type"] = type2string(system)
+    system_data["particle_spacing"] = particle_spacing(system, 1)
+    system_data["adhesion_coefficient"] = system.adhesion_coefficient
+    add_system_data!(system_data, system.boundary_model)
+    add_system_data!(system_data, system.movement)
+end
+
+function add_system_data!(system_data, system::BoundaryDEMSystem)
+    system_data["system_type"] = type2string(system)
+    system_data["particle_spacing"] = particle_spacing(system, 1)
+    system_data["normal_stiffness"] = system.normal_stiffness
+end
+
+function add_system_data!(system_data, system::DEMSystem)
+    system_data["system_type"] = type2string(system)
+    system_data["particle_spacing"] = particle_spacing(system, 1)
+    system_data["damping_coefficient"] = system.damping_coefficient
+    system_data["acceleration"] = system.acceleration
+    add_system_data!(system_data, system.contact_model)
+end
+
+function add_system_data!(system_data, system::OpenBoundarySPHSystem)
+    system_data["system_type"] = type2string(system)
+    system_data["fluid_system_index"] = system.fluid_system_index[]
+    system_data["smoothing_length"] = system.smoothing_length
+    add_system_data!(system_data, system.boundary_model)
+
+    system_data["boundary_zones"] = Dict{String, Any}()
+    for (indice, boundary_zone) in enumerate(system.boundary_zones)
+        add_system_data!(system_data["boundary_zones"], boundary_zone, indice)
+    end
+end
+
+function add_system_data!(system_data, system::ParticlePackingSystem)
+    system_data["system_type"] = type2string(system)
+    system_data["particle_spacing"] = system.particle_spacing
+    system_data["smoothing_kernel"] = type2string(system.smoothing_kernel)
+    system_data["smoothing_length_interpolation"] = system.smoothing_length_interpolation
+    system_data["background_pressure"] = system.background_pressure
+    system_data["place_on_shell"] = system.place_on_shell
+    system_data["shift_length"] = system.shift_length
+end
+
+function add_system_data!(system_data, boundary_model::BoundaryModelDummyParticles)
+    system_data["boundary_model"] = Dict{String, Any}()
+    system_data["boundary_model"]["model"] = type2string(boundary_model)
+    system_data["boundary_model"]["smoothing_kernel"] = type2string(boundary_model.smoothing_kernel)
+    system_data["boundary_model"]["smoothing_length"] = boundary_model.smoothing_length
+    system_data["boundary_model"]["density_calculator"] = type2string(boundary_model.density_calculator)
+    add_system_data!(system_data["boundary_model"], boundary_model.state_equation)
+    add_system_data!(system_data["boundary_model"], boundary_model.viscosity)
+    add_system_data!(system_data["boundary_model"], boundary_model.correction)
+end
+
+function add_system_data!(system_data, boundary_model::BoundaryModelMonaghanKajtar)
+    system_data["boundary_model"] = Dict{String, Any}()
+    system_data["boundary_model"]["model"] = type2string(boundary_model)
+    system_data["boundary_model"]["beta"] = boundary_model.beta
+    system_data["boundary_model"]["K"] = boundary_model.K
+    add_system_data!(system_data["boundary_model"], boundary_model.viscosity)
+end
+
+function add_system_data!(system_data, boundary_model::BoundaryModelMirroringTafuni)
+    system_data["boundary_model"] = Dict{String, Any}()
+    system_data["boundary_model"]["model"] = type2string(boundary_model)
+    system_data["boundary_model"]["mirror_method"] = type2string(boundary_model.mirror_method)
+end
+
+function add_system_data!(system_data, boundary_model::BoundaryModelCharacteristicsLastiwka)
+    system_data["boundary_model"] = Dict{String, Any}()
+    system_data["boundary_model"]["model"] = type2string(boundary_model)
+    system_data["boundary_model"]["extrapolate_reference_values"] = boundary_model.extrapolate_reference_values
+end
+
+function add_system_data!(system_data, contact_model::HertzContactModel)
+    system_data["contact_model"] = Dict{String, Any}()
+    system_data["contact_model"]["model"] = type2string(contact_model)
+    system_data["contact_model"]["elastic_modulus"] = contact_model.elastic_modulus
+    system_data["contact_model"]["poissons_ratio"] = contact_model.poissons_ratio
+end
+
+function add_system_data!(system_data, contact_model::LinearContactModel)
+    system_data["contact_model"] = Dict{String, Any}()
+    system_data["contact_model"]["model"] = type2string(contact_model)
+    system_data["contact_model"]["normal_stiffness"] = contact_model.normal_stiffness
+end
+
+function add_system_data!(system_data, state_equation::StateEquationCole)
+    system_data["state_equation"] = Dict{String, Any}()
+    system_data["state_equation"]["model"] = type2string(state_equation)
+    system_data["state_equation"]["reference_density"] = state_equation.reference_density
+    system_data["state_equation"]["background_pressure"] = state_equation.background_pressure
+    system_data["state_equation"]["exponent"] = state_equation.exponent
+end
+
+function add_system_data!(system_data, state_equation::StateEquationIdealGas)
+    system_data["state_equation"] = Dict{String, Any}()
+    system_data["state_equation"]["model"] = type2string(state_equation)
+    system_data["state_equation"]["reference_density"] = state_equation.reference_density
+    system_data["state_equation"]["background_pressure"] = state_equation.background_pressure
+    system_data["state_equation"]["gamma"] = state_equation.gamma
+end
+
+function add_system_data!(system_data, viscosity::Union{ViscosityAdami, ViscosityMorris})
+    system_data["viscosity_model"] = Dict{String, Any}()
+    system_data["viscosity_model"]["model"] = type2string(viscosity)
+    system_data["viscosity_model"]["nu"] = viscosity.nu
+    system_data["viscosity_model"]["epsilon"] = viscosity.epsilon
+end
+
+function add_system_data!(system_data,
+                          viscosity::Union{ViscosityAdamiSGS, ViscosityMorrisSGS})
+    system_data["viscosity_model"] = Dict{String, Any}()
+    system_data["viscosity_model"]["model"] = type2string(viscosity)
+    system_data["viscosity_model"]["nu"] = viscosity.nu
+    system_data["viscosity_model"]["C_S"] = viscosity.C_S
+    system_data["viscosity_model"]["epsilon"] = viscosity.epsilon
+end
+
+function add_system_data!(system_data, viscosity::ArtificialViscosityMonaghan)
+    system_data["viscosity_model"] = Dict{String, Any}()
+    system_data["viscosity_model"]["model"] = type2string(viscosity)
+    system_data["viscosity_model"]["alpha"] = viscosity.alpha
+    system_data["viscosity_model"]["beta"] = viscosity.beta
+    system_data["viscosity_model"]["epsilon"] = viscosity.epsilon
+end
+
+function add_system_data!(system_data,
+                          density_diffusion::Union{DensityDiffusionAntuono,
+                                                   DensityDiffusionMolteniColagrossi})
+    system_data["density_diffusion"] = Dict{String, Any}()
+    system_data["density_diffusion"]["model"] = type2string(density_diffusion)
+    system_data["density_diffusion"]["delta"] = density_diffusion.delta
+end
+
+function add_system_data!(system_data, density_diffusion::DensityDiffusionFerrari)
+    system_data["density_diffusion"] = Dict{String, Any}()
+    system_data["density_diffusion"]["model"] = type2string(density_diffusion)
+end
+
+function add_system_data!(system_data, correction::AkinciFreeSurfaceCorrection)
+    system_data["correction_method"] = Dict{String, Any}()
+    system_data["correction_method"]["model"] = type2string(correction)
+    system_data["correction_method"]["rho0"] = correction.rho0
+end
+
+function add_system_data!(system_data,
+                          correction::Union{BlendedGradientCorrection, GradientCorrection,
+                                            KernelCorrection, MixedKernelGradientCorrection,
+                                            ShepardKernelCorrection})
+    system_data["correction_method"] = Dict{String, Any}()
+    system_data["correction_method"]["model"] = type2string(correction)
+end
+
+function add_system_data!(system_data,
+                          surface_tension::Union{CohesionForceAkinci, SurfaceTensionAkinci,
+                                                 SurfaceTensionMorris,
+                                                 SurfaceTensionMomentumMorris})
+    system_data["surface_tension"] = Dict{String, Any}()
+    system_data["surface_tension"]["model"] = type2string(surface_tension)
+    system_data["surface_tension"]["surface_tension_coefficient"] = surface_tension.surface_tension_coefficient
+end
+
+function add_system_data!(system_data, surface_normal_method::ColorfieldSurfaceNormal)
+    system_data["surface_normal_method"] = Dict{String, Any}()
+    system_data["surface_normal_method"]["model"] = type2string(surface_normal_method)
+    system_data["surface_normal_method"]["boundary_contact_threshold"] = surface_normal_method.boundary_contact_threshold
+    system_data["surface_normal_method"]["ideal_density_threshold"] = surface_normal_method.ideal_density_threshold
+end
+
+function add_system_data!(system_data, boundary_zone::BoundaryZone, indice)
+    zone_name = "boundary_zone_" * string(indice)
+    system_data[zone_name] = Dict{String, Any}()
+    system_data[zone_name]["spanning_set"] = boundary_zone.spanning_set
+    system_data[zone_name]["zone_origin"] = boundary_zone.zone_origin
+    system_data[zone_name]["zone_width"] = boundary_zone.zone_width
+    system_data[zone_name]["flow_direction"] = boundary_zone.flow_direction
+    system_data[zone_name]["plane_normal"] = boundary_zone.plane_normal
+    system_data[zone_name]["reference_values"] = boundary_zone.reference_values
+    system_data[zone_name]["average_inflow_velocity"] = boundary_zone.average_inflow_velocity
+    system_data[zone_name]["prescribed_density"] = boundary_zone.prescribed_density
+    system_data[zone_name]["prescribed_pressure"] = boundary_zone.prescribed_pressure
+    system_data[zone_name]["prescribed_velocity"] = boundary_zone.prescribed_velocity
+end
+
+function add_system_data!(system_data, movement::BoundaryMovement)
+    system_data["movement"] = Dict{String, Any}()
+    system_data["movement"]["model"] = type2string(movement)
+    system_data["movement"]["movement_function"] = type2string(movement.movement_function)
+    system_data["movement"]["is_moving"] = type2string(movement.is_moving)
+    system_data["movement"]["moving_particles"] = movement.moving_particles
+end
+
+function add_system_data!(system_data, penalty_force::PenaltyForceGanzenmueller)
+    system_data["penalty_force"] = Dict{String, Any}()
+    system_data["penalty_force"]["model"] = type2string(penalty_force)
+    system_data["penalty_force"]["alpha"] = penalty_force.alpha
+end
+
+function add_system_data!(system_data, shifting_technique::TransportVelocityAdami)
+    system_data["shifting_technique"] = Dict{String, Any}()
+    system_data["shifting_technique"]["model"] = type2string(shifting_technique)
+    system_data["shifting_technique"]["background_pressure"] = shifting_technique.background_pressure
+end
+
+function add_system_data!(system_data, shifting_technique::ParticleShiftingTechnique)
+    system_data["shifting_technique"] = Dict{String, Any}()
+    system_data["shifting_technique"]["model"] = type2string(shifting_technique)
+end
diff --git a/src/io/write_vtk.jl b/src/io/write_vtk.jl
index c17410d698..6f0b772189 100644
--- a/src/io/write_vtk.jl
+++ b/src/io/write_vtk.jl
@@ -10,7 +10,7 @@ end
 
 """
     trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out", prefix="",
-              write_meta_data=true, max_coordinates=Inf, custom_quantities...)
+              max_coordinates=Inf, custom_quantities...)
 
 Convert Trixi simulation data to VTK format.
 
@@ -25,7 +25,6 @@ Convert Trixi simulation data to VTK format.
                             separate files. This number is just appended to the filename.
 - `output_directory="out"`: Output directory path.
 - `prefix=""`:              Prefix for output files.
-- `write_meta_data=true`:   Write meta data.
 - `max_coordinates=Inf`     The coordinates of particles will be clipped if their absolute
                             values exceed this threshold.
 - `custom_quantities...`:   Additional custom quantities to include in the VTK output.
@@ -49,20 +48,19 @@ trixi2vtk(sol.u[end], semi, 0.0, iter=1, my_custom_quantity=kinetic_energy)
 ```
 """
 function trixi2vtk(vu_ode, semi, t; iter=nothing, output_directory="out",
-                   prefix="", write_meta_data=true, git_hash=compute_git_hash(),
-                   max_coordinates=Inf, custom_quantities...)
+                   prefix="", git_hash=compute_git_hash(), max_coordinates=Inf,
+                   custom_quantities...)
 
     # The first argument is not necessary in most cases. Since it is usually not available to the user,
     # this API wrapper makes it optional.
     # Note that custom quantities using the fluid acceleration will not work and return NaN acceleration.
     return trixi2vtk(fill!(similar(vu_ode), NaN), vu_ode, semi, t; iter, output_directory,
-                     prefix, write_meta_data, git_hash, max_coordinates,
-                     custom_quantities...)
+                     prefix, git_hash, max_coordinates, custom_quantities...)
 end
 
 function trixi2vtk(dvdu_ode, vu_ode, semi, t; iter=nothing, output_directory="out",
-                   prefix="", write_meta_data=true, git_hash=compute_git_hash(),
-                   max_coordinates=Inf, custom_quantities...)
+                   prefix="", git_hash=compute_git_hash(), max_coordinates=Inf,
+                   custom_quantities...)
     (; systems) = semi
 
     # Update quantities that are stored in the systems. These quantities (e.g. pressure)
@@ -81,14 +79,14 @@ function trixi2vtk(dvdu_ode, vu_ode, semi, t; iter=nothing, output_directory="ou
 
         trixi2vtk(system, dvdu_ode, vu_ode, semi, t, periodic_box;
                   system_name=filenames[system_index], output_directory, iter, prefix,
-                  write_meta_data, git_hash, max_coordinates, custom_quantities...)
+                  git_hash, max_coordinates, custom_quantities...)
     end
 end
 
 # Convert data for a single TrixiParticle system to VTK format
 function trixi2vtk(system_, dvdu_ode_, vu_ode_, semi_, t, periodic_box;
                    output_directory="out", prefix="", iter=nothing,
-                   system_name=vtkname(system_), write_meta_data=true, max_coordinates=Inf,
+                   system_name=vtkname(system_), max_coordinates=Inf,
                    git_hash=compute_git_hash(), custom_quantities...)
     mkpath(output_directory)
 
@@ -105,16 +103,12 @@ function trixi2vtk(system_, dvdu_ode_, vu_ode_, semi_, t, periodic_box;
     v = wrap_v(v_ode, system, semi)
     u = wrap_u(u_ode, system, semi)
 
-    # handle "_" on optional pre/postfix strings
-    add_opt_str_pre(str) = (str === "" ? "" : "$(str)_")
-    add_opt_str_post(str) = (str === nothing ? "" : "_$(str)")
-
     file = joinpath(output_directory,
-                    add_opt_str_pre(prefix) * "$system_name"
-                    * add_opt_str_post(iter))
+                    add_underscore_to_optional_prefix(prefix) * "$system_name"
+                    * add_underscore_to_optional_postfix(iter))
 
     collection_file = joinpath(output_directory,
-                               add_opt_str_pre(prefix) * "$system_name")
+                               add_underscore_to_optional_prefix(prefix) * "$system_name")
 
     # Reset the collection when the iteration is 0
     pvd = paraview_collection(collection_file; append=iter > 0)
@@ -134,7 +128,7 @@ function trixi2vtk(system_, dvdu_ode_, vu_ode_, semi_, t, periodic_box;
 
     @trixi_timeit timer() "write to vtk" vtk_grid(file, points, cells) do vtk
         # Dispatches based on the different system types e.g. FluidSystem, TotalLagrangianSPHSystem
-        write2vtk!(vtk, v, u, t, system, write_meta_data=write_meta_data)
+        write2vtk!(vtk, v, u, t, system)
 
         # Store particle index
         vtk["index"] = eachparticle(system)
@@ -142,12 +136,7 @@ function trixi2vtk(system_, dvdu_ode_, vu_ode_, semi_, t, periodic_box;
         vtk["ndims"] = ndims(system)
 
         vtk["particle_spacing"] = [particle_spacing(system, particle)
-                                   for particle in eachparticle(system)]
-
-        if write_meta_data
-            vtk["solver_version"] = git_hash
-            vtk["julia_version"] = string(VERSION)
-        end
+                                   for particle in active_particles(system)]
 
         # Extract custom quantities for this system
         if !isempty(custom_quantities)
@@ -285,13 +274,13 @@ function trixi2vtk(initial_condition::InitialCondition; output_directory="out",
                      pressure=pressure, custom_quantities...)
 end
 
-function write2vtk!(vtk, v, u, t, system; write_meta_data=true)
+function write2vtk!(vtk, v, u, t, system)
     vtk["velocity"] = view(v, 1:ndims(system), :)
 
     return vtk
 end
 
-function write2vtk!(vtk, v, u, t, system::DEMSystem; write_meta_data=true)
+function write2vtk!(vtk, v, u, t, system::DEMSystem)
     vtk["velocity"] = view(v, 1:ndims(system), :)
     vtk["mass"] = [hydrodynamic_mass(system, particle)
                    for particle in eachparticle(system)]
@@ -300,7 +289,7 @@ function write2vtk!(vtk, v, u, t, system::DEMSystem; write_meta_data=true)
     return vtk
 end
 
-function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true)
+function write2vtk!(vtk, v, u, t, system::FluidSystem)
     vtk["velocity"] = [current_velocity(v, system, particle)
                        for particle in eachparticle(system)]
     vtk["density"] = [current_density(v, system, particle)
@@ -334,9 +323,14 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true)
             surface_tension[1:ndims(system),
                             particle] .+= surface_tension_force(surface_tension_a,
                                                                 surface_tension_b,
-                                                                system, system, particle,
-                                                                neighbor, pos_diff,
-                                                                distance, rho_a, rho_b,
+                                                                system,
+                                                                system,
+                                                                particle,
+                                                                neighbor,
+                                                                pos_diff,
+                                                                distance,
+                                                                rho_a,
+                                                                rho_b,
                                                                 grad_kernel)
         end
         vtk["surface_tension"] = surface_tension
@@ -349,41 +343,6 @@ function write2vtk!(vtk, v, u, t, system::FluidSystem; write_meta_data=true)
         end
     end
 
-    if write_meta_data
-        vtk["acceleration"] = system.acceleration
-        vtk["viscosity"] = type2string(system.viscosity)
-        write2vtk!(vtk, system.viscosity)
-        vtk["smoothing_kernel"] = type2string(system.smoothing_kernel)
-        vtk["smoothing_length_factor"] = system.cache.smoothing_length_factor
-        vtk["density_calculator"] = type2string(system.density_calculator)
-
-        if system isa WeaklyCompressibleSPHSystem
-            vtk["solver"] = "WCSPH"
-
-            vtk["correction_method"] = type2string(system.correction)
-            if system.correction isa AkinciFreeSurfaceCorrection
-                vtk["correction_rho0"] = system.correction.rho0
-            end
-
-            if system.state_equation isa StateEquationCole
-                vtk["state_equation_exponent"] = system.state_equation.exponent
-            end
-
-            if system.state_equation isa StateEquationIdealGas
-                vtk["state_equation_gamma"] = system.state_equation.gamma
-            end
-
-            vtk["state_equation"] = type2string(system.state_equation)
-            vtk["state_equation_rho0"] = system.state_equation.reference_density
-            vtk["state_equation_pa"] = system.state_equation.background_pressure
-            vtk["state_equation_c"] = system.state_equation.sound_speed
-            vtk["solver"] = "WCSPH"
-        else
-            vtk["solver"] = "EDAC"
-            vtk["sound_speed"] = system.sound_speed
-        end
-    end
-
     return vtk
 end
 
@@ -402,7 +361,7 @@ function write2vtk!(vtk, viscosity::ArtificialViscosityMonaghan)
     vtk["viscosity_epsilon"] = viscosity.epsilon
 end
 
-function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_data=true)
+function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem)
     n_fixed_particles = nparticles(system) - n_moving_particles(system)
 
     vtk["velocity"] = [current_velocity(v, system, particle)
@@ -416,6 +375,11 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d
                            initial_coords(system, particle)
                            for particle in eachparticle(system)]
 
+    vtk["lame_lambda"] = system.lame_lambda
+    vtk["lame_mu"] = system.lame_mu
+    vtk["young_modulus"] = system.young_modulus
+    vtk["poisson_ratio"] = system.poisson_ratio
+
     sigma = cauchy_stress(system)
     vtk["sigma_11"] = sigma[1, 1, :]
     vtk["sigma_22"] = sigma[2, 2, :]
@@ -425,18 +389,10 @@ function write2vtk!(vtk, v, u, t, system::TotalLagrangianSPHSystem; write_meta_d
 
     vtk["material_density"] = system.material_density
 
-    if write_meta_data
-        vtk["lame_lambda"] = system.lame_lambda
-        vtk["lame_mu"] = system.lame_mu
-        vtk["smoothing_kernel"] = type2string(system.smoothing_kernel)
-        vtk["smoothing_length_factor"] = initial_smoothing_length(system) /
-                                         particle_spacing(system, 1)
-    end
-
-    write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data)
+    write2vtk!(vtk, v, u, t, system.boundary_model, system)
 end
 
-function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data=true)
+function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem)
     vtk["velocity"] = [current_velocity(v, system, particle)
                        for particle in eachparticle(system)]
     vtk["density"] = [current_density(v, system, particle)
@@ -447,34 +403,19 @@ function write2vtk!(vtk, v, u, t, system::OpenBoundarySPHSystem; write_meta_data
     return vtk
 end
 
-function write2vtk!(vtk, v, u, t, system::BoundarySPHSystem; write_meta_data=true)
-    write2vtk!(vtk, v, u, t, system.boundary_model, system, write_meta_data=write_meta_data)
+function write2vtk!(vtk, v, u, t, system::BoundarySPHSystem)
+    write2vtk!(vtk, v, u, t, system.boundary_model, system)
 end
 
-function write2vtk!(vtk, v, u, t, model::Nothing, system; write_meta_data=true)
+function write2vtk!(vtk, v, u, t, model::Nothing, system)
     return vtk
 end
 
-function write2vtk!(vtk, v, u, t, model::BoundaryModelMonaghanKajtar, system;
-                    write_meta_data=true)
-    if write_meta_data
-        vtk["boundary_model"] = "BoundaryModelMonaghanKajtar"
-        vtk["boundary_spacing_ratio"] = model.beta
-        vtk["boundary_K"] = model.K
-    end
+function write2vtk!(vtk, v, u, t, model::BoundaryModelMonaghanKajtar, system)
+    return vtk
 end
 
-function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system;
-                    write_meta_data=true)
-    if write_meta_data
-        vtk["boundary_model"] = "BoundaryModelDummyParticles"
-        vtk["smoothing_kernel"] = type2string(model.smoothing_kernel)
-        vtk["smoothing_length"] = model.smoothing_length
-        vtk["density_calculator"] = type2string(model.density_calculator)
-        vtk["state_equation"] = type2string(model.state_equation)
-        vtk["viscosity_model"] = type2string(model.viscosity)
-    end
-
+function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system)
     vtk["hydrodynamic_density"] = current_density(v, system)
     vtk["pressure"] = model.pressure
 
@@ -489,6 +430,6 @@ function write2vtk!(vtk, v, u, t, model::BoundaryModelDummyParticles, system;
     end
 end
 
-function write2vtk!(vtk, v, u, t, system::BoundaryDEMSystem; write_meta_data=true)
+function write2vtk!(vtk, v, u, t, system::BoundaryDEMSystem)
     return vtk
 end
diff --git a/src/preprocessing/particle_packing/system.jl b/src/preprocessing/particle_packing/system.jl
index 100d0f4963..42502e9afc 100644
--- a/src/preprocessing/particle_packing/system.jl
+++ b/src/preprocessing/particle_packing/system.jl
@@ -216,12 +216,9 @@ end
 
 @inline requires_update_callback(system::ParticlePackingSystem) = true
 
-function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem; write_meta_data=true)
+function write2vtk!(vtk, v, u, t, system::ParticlePackingSystem)
     vtk["velocity"] = [advection_velocity(v, system, particle)
                        for particle in eachparticle(system)]
-    if write_meta_data
-        vtk["signed_distances"] = system.signed_distances
-    end
 end
 
 # Skip for fixed systems