diff --git a/SPlisHSPlasH/CMakeLists.txt b/SPlisHSPlasH/CMakeLists.txt
index c47f1f1d..a11d1d22 100644
--- a/SPlisHSPlasH/CMakeLists.txt
+++ b/SPlisHSPlasH/CMakeLists.txt
@@ -195,6 +195,8 @@ add_library(SPlisHSPlasH
 	AnimationField.h
 	AnimationFieldSystem.h
 	AnimationFieldSystem.cpp
+	DynamicParameterSystem.h
+	DynamicParameterSystem.cpp
 	BoundaryModel.cpp
 	BoundaryModel.h
 	BoundaryModel_Akinci2012.cpp
diff --git a/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp b/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp
index 9ab87b34..f9046447 100644
--- a/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp
+++ b/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp
@@ -237,10 +237,11 @@ void TimeStepDFSPH::step()
 	}
 
 	//////////////////////////////////////////////////////////////////////////
-	// emit new particles and perform an animation field step
+	// emit new particles, perform an animation field step and update parameters
 	//////////////////////////////////////////////////////////////////////////
 	sim->emitParticles();
 	sim->animateParticles();
+	sim->updateDynamicParameters();
 
 	//////////////////////////////////////////////////////////////////////////
 	// Compute new time
diff --git a/SPlisHSPlasH/DynamicParameterSystem.cpp b/SPlisHSPlasH/DynamicParameterSystem.cpp
new file mode 100644
index 00000000..f420fb56
--- /dev/null
+++ b/SPlisHSPlasH/DynamicParameterSystem.cpp
@@ -0,0 +1,410 @@
+#include "DynamicParameterSystem.h"
+#include "FluidModel.h"
+#include "Simulation.h"
+#include "TimeManager.h"
+#include "Utilities/Logger.h"
+#include "WCSPH/TimeStepWCSPH.h"
+#include "PF/TimeStepPF.h"
+#include "ICSPH/TimeStepICSPH.h"
+#include "Viscosity/Viscosity_Standard.h"
+#include "Viscosity/Viscosity_Bender2017.h"
+#include "Viscosity/Viscosity_Peer2015.h"
+#include "Viscosity/Viscosity_Peer2016.h"
+#include "Viscosity/Viscosity_Takahashi2015.h"
+#include "Viscosity/Viscosity_Weiler2018.h"
+#include "SurfaceTension/SurfaceTension_Becker2007.h"
+#include "SurfaceTension/SurfaceTension_Akinci2013.h"
+#include "SurfaceTension/SurfaceTension_He2014.h"
+#include "SurfaceTension/SurfaceTension_Jeske2023.h"
+#include "SurfaceTension/SurfaceTension_ZorillaRitter2020.h"
+#include <algorithm>
+
+using namespace SPH;
+
+DynamicParameterSystem::DynamicParameterSystem() :
+    m_t_double(0.0),
+    m_dt_double(0.0)
+{
+    LOG_INFO << "Dynamic Parameters created"; //Debug
+}
+
+DynamicParameterSystem::~DynamicParameterSystem()
+{
+    reset();
+}
+
+void DynamicParameterSystem::step()
+{
+    const Real currentTime = TimeManager::getCurrent()->getTime();
+    const Real dt = TimeManager::getCurrent()->getTimeStepSize();
+
+    m_t_double = currentTime;
+    m_dt_double = dt;
+
+    for (auto& schedule : m_schedules)
+    {
+        if (!schedule.active)
+            continue;
+
+        Real new_value;
+
+        if (schedule.use_expression && schedule.compiled_expr.get() != nullptr)
+        {
+            // Evaluate
+            new_value = static_cast<Real>(te_eval(schedule.compiled_expr.get()));
+            if (std::isnan(new_value) || std::isinf(new_value)) {
+                LOG_ERR << "Expression evaluation failed for " << schedule.parameter_name;
+                schedule.active = false;
+                continue;
+            }
+        }
+        else
+        {
+            if (schedule.use_step_function) 
+            {
+                new_value = stepTimeline(schedule, currentTime);
+                if (std::isnan(new_value))
+                    continue;
+            }  
+            else
+            {
+                new_value = interpolateTimeline(schedule, currentTime);
+                if (std::isnan(new_value))
+                    continue;
+            }
+            if (currentTime >= schedule.timeline.back().first)
+                schedule.active = false;
+        }
+
+        schedule.current_value = new_value;
+        applyParameterValue(schedule.fluidId, schedule.parameter_name, new_value);
+    }
+}
+
+void DynamicParameterSystem::reset()
+{
+    for (auto& schedule : m_schedules)
+    {
+        applyParameterValue(schedule.fluidId, schedule.parameter_name, schedule.default_value);
+    }
+    m_schedules.clear();
+}
+
+bool DynamicParameterSystem::addTimelineSchedule(const std::string& fluidId, const std::string& paramName, const std::vector<std::pair<Real, Real>>& timeline, Real defaultValue, bool stepFunction)
+{
+    ParameterSchedule schedule;
+    schedule.fluidId = fluidId;
+    schedule.parameter_name = paramName;
+    schedule.timeline = timeline;
+    schedule.default_value = defaultValue;
+    schedule.current_value = defaultValue;
+    schedule.use_expression = false;
+    schedule.use_step_function = stepFunction;
+    schedule.active = true;
+
+    // Order timeline
+    std::sort(schedule.timeline.begin(), schedule.timeline.end());
+
+    m_schedules.push_back(std::move(schedule));
+
+    LOG_INFO << "Added timeline schedule for " << paramName << " in Fluid " << fluidId
+        << " with " << timeline.size() << " points"
+        << (stepFunction ? " (step function)" : " (interpolated)");
+    return true;
+}
+
+bool DynamicParameterSystem::addExpressionSchedule(const std::string& fluidId, const std::string& paramName, const std::string& expression, Real defaultValue)
+{
+    ParameterSchedule schedule;
+    schedule.fluidId = fluidId;
+    schedule.parameter_name = paramName;
+    schedule.expression = expression;
+    schedule.default_value = defaultValue;
+    schedule.current_value = defaultValue;
+    schedule.use_expression = true;
+    schedule.active = true;
+
+    if (!compileExpression(schedule))
+    {
+        LOG_ERR << "Failed to compile expression for " << paramName
+            << ": " << expression;
+        return false;
+    }
+
+    m_schedules.push_back(std::move(schedule)); // why the move here?
+
+    LOG_INFO << "Added expression schedule for " << paramName << " in Fluid " << fluidId << ": " << expression;
+    return true;
+}
+
+bool DynamicParameterSystem::compileExpression(ParameterSchedule& schedule)
+{
+    te_variable vars[] = {
+        {"t", &m_t_double},
+        {"dt", &m_dt_double}
+    };
+    const int numVars = 2;
+
+    int err;
+    te_expr* raw_expr = te_compile(schedule.expression.c_str(), vars, numVars, &err);
+
+    if (raw_expr != nullptr && err == 0) {
+        schedule.compiled_expr.reset(raw_expr);
+    }
+
+    return (schedule.compiled_expr != nullptr && err == 0);
+}
+
+Real DynamicParameterSystem::interpolateTimeline(ParameterSchedule& schedule, Real currentTime)
+{
+    const auto& timeline = schedule.timeline;
+
+    if (timeline.empty())
+        return schedule.default_value;
+
+    // Before first
+    if (currentTime < timeline.front().first)
+    {
+        schedule.current_timeline_index = -1;
+        return std::numeric_limits<Real>::quiet_NaN();
+    }
+
+    // After last
+    if (currentTime >= timeline.back().first)
+    {
+        return timeline.back().second;
+    }
+
+    int old_index = schedule.current_timeline_index;
+
+    while (schedule.current_timeline_index < static_cast<int>(timeline.size()) - 1 && currentTime > timeline[schedule.current_timeline_index + 1].first)
+    {
+        schedule.current_timeline_index++;
+    }
+
+    // Logging
+    if (schedule.current_timeline_index != old_index)
+    {
+        int i = schedule.current_timeline_index;
+        LOG_INFO << schedule.parameter_name << " (Fluid: " << schedule.fluidId
+            << ") will interpolate from " << timeline[i].second
+            << " to " << timeline[i + 1].second
+            << " between t=" << timeline[i].first
+            << " and t=" << timeline[i + 1].first;
+    }
+
+    // Interpolate
+    int i = schedule.current_timeline_index;
+    const Real t1 = timeline[i].first;
+    const Real v1 = timeline[i].second;
+    const Real t2 = timeline[i + 1].first;
+    const Real v2 = timeline[i + 1].second;
+
+    const Real alpha = (currentTime - t1) / (t2 - t1);
+    return v1 + alpha * (v2 - v1);
+}
+
+Real DynamicParameterSystem::stepTimeline(ParameterSchedule& schedule, Real currentTime)
+{
+    const auto& timeline = schedule.timeline;
+    if (timeline.empty())
+        return schedule.default_value;
+
+    // Before first
+    if (currentTime < timeline.front().first)
+    {
+        schedule.current_timeline_index = -1;
+        return std::numeric_limits<Real>::quiet_NaN();
+    }
+
+    int old_index = schedule.current_timeline_index;
+
+    while (schedule.current_timeline_index < static_cast<int>(timeline.size()) - 1 && currentTime >= timeline[schedule.current_timeline_index + 1].first)
+    {
+        schedule.current_timeline_index++;
+    }
+
+    // Logging
+    if (schedule.current_timeline_index != old_index)
+    {
+        int i = schedule.current_timeline_index;
+        LOG_INFO << schedule.parameter_name << " (Fluid: " << schedule.fluidId
+            << ") changed to " << timeline[i].second
+            << " at t=" << timeline[i].first;
+    }
+
+    return timeline[schedule.current_timeline_index].second;
+}
+
+void DynamicParameterSystem::applyParameterValue(const std::string& fluidId, const std::string& paramName, Real value)
+{
+    const Real currentTime = TimeManager::getCurrent()->getTime();
+    Simulation* sim = Simulation::getCurrent();
+
+    if (paramName == "stiffness")
+    {
+        if (sim->getSimulationMethod() == 0) // WCSPH
+        {
+            TimeStepWCSPH* timeStep = static_cast<TimeStepWCSPH*>(sim->getTimeStep());
+            timeStep->setValue(TimeStepWCSPH::STIFFNESS, value);
+            return;
+        } 
+        else if (sim->getSimulationMethod() == 5) // PF
+        {
+            TimeStepPF* timeStep = static_cast<TimeStepPF*>(sim->getTimeStep());
+            timeStep->setValue(TimeStepPF::STIFFNESS, value);
+            return;
+        } else {
+            LOG_WARN << "Not a valid simulation method.";
+        }
+    }
+
+    if (paramName == "lambda")
+    {
+        if (sim->getSimulationMethod() == 6) // ICSPH
+        {
+            TimeStepICSPH* timeStep = static_cast<TimeStepICSPH*>(sim->getTimeStep());
+            timeStep->setValue(TimeStepICSPH::LAMBDA, value);
+            return;
+        }
+        else {
+            LOG_WARN << "Not a valid simulation method.";
+        }
+    }
+
+    for (unsigned int i = 0; i < sim->numberOfFluidModels(); i++)
+    {
+        if (sim->getFluidModel(i)->getId() == fluidId)
+        {
+            FluidModel* model = sim->getFluidModel(i);
+            if (paramName == "density0")
+            {
+                model->setDensity0(value);
+            }
+            else if (paramName == "viscosity")
+            {
+                if (model->getViscosityMethod() == 1) // Standard
+                {
+                    Viscosity_Standard* viscMethod = static_cast<Viscosity_Standard*>(model->getViscosityBase());
+                    viscMethod->setValue(Viscosity_Standard::VISCOSITY_COEFFICIENT, value);
+                }
+                else if (model->getViscosityMethod() == 2) // Bender 2017
+                {
+                    Viscosity_Bender2017* viscMethod = static_cast<Viscosity_Bender2017*>(model->getViscosityBase());
+                    viscMethod->setValue(Viscosity_Bender2017::VISCOSITY_COEFFICIENT, value);
+                }
+                else if (model->getViscosityMethod() == 3) // Peer 2015
+                {
+                    Viscosity_Peer2015* viscMethod = static_cast<Viscosity_Peer2015*>(model->getViscosityBase());
+                    viscMethod->setValue(Viscosity_Peer2015::VISCOSITY_COEFFICIENT, value);
+                }
+                else if (model->getViscosityMethod() == 4) // Peer 2016
+                {
+                    Viscosity_Peer2016* viscMethod = static_cast<Viscosity_Peer2016*>(model->getViscosityBase());
+                    viscMethod->setValue(Viscosity_Peer2016::VISCOSITY_COEFFICIENT, value);
+                }
+                else if (model->getViscosityMethod() == 5) // Takahashi 2015
+                {
+                    Viscosity_Takahashi2015* viscMethod = static_cast<Viscosity_Takahashi2015*>(model->getViscosityBase());
+                    viscMethod->setValue(Viscosity_Takahashi2015::VISCOSITY_COEFFICIENT, value);
+                }
+                else if (model->getViscosityMethod() == 6) // Weiler 2018
+                {
+                    Viscosity_Weiler2018* viscMethod = static_cast<Viscosity_Weiler2018*>(model->getViscosityBase());
+                    viscMethod->setValue(Viscosity_Weiler2018::VISCOSITY_COEFFICIENT, value);
+                }
+            }
+            else if (paramName == "viscosityBoundary")
+            {
+                if (model->getViscosityMethod() == 1) // Standard
+                {
+                    Viscosity_Standard* viscMethod = static_cast<Viscosity_Standard*>(model->getViscosityBase());
+                    viscMethod->setValue(Viscosity_Standard::VISCOSITY_COEFFICIENT_BOUNDARY, value);
+                }
+                else if (model->getViscosityMethod() == 6) // Weiler 2018
+                {
+                    Viscosity_Weiler2018* viscMethod = static_cast<Viscosity_Weiler2018*>(model->getViscosityBase());
+                    viscMethod->setValue(Viscosity_Weiler2018::VISCOSITY_COEFFICIENT_BOUNDARY, value);
+                }
+            }
+            else if (paramName == "surfaceTension")
+            {
+                if (model->getSurfaceTensionMethod() == 1) // Becker & Teschner 2007
+                {
+                    SurfaceTension_Becker2007* stMethod = static_cast<SurfaceTension_Becker2007*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_Becker2007::SURFACE_TENSION, value);
+                }
+                else if (model->getSurfaceTensionMethod() == 2) // Akinci 2013
+                {
+                    SurfaceTension_Akinci2013* stMethod = static_cast<SurfaceTension_Akinci2013*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_Akinci2013::SURFACE_TENSION, value);
+                }
+                else if (model->getSurfaceTensionMethod() == 3) // He 2014
+                {
+                    SurfaceTension_He2014* stMethod = static_cast<SurfaceTension_He2014*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_He2014::SURFACE_TENSION, value);
+                }
+                else if (model->getSurfaceTensionMethod() == 4) // Jeske 2023
+                {
+                    SurfaceTension_Jeske2023* stMethod = static_cast<SurfaceTension_Jeske2023*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_Jeske2023::SURFACE_TENSION, value);
+                }
+                else {
+                    LOG_WARN << "Not a valid surface tension method.";
+                }
+            }
+            else if (paramName == "surfaceTensionBoundary")
+            {
+                if (model->getSurfaceTensionMethod() == 1) // Becker & Teschner 2007
+                {
+                    SurfaceTension_Becker2007* stMethod = static_cast<SurfaceTension_Becker2007*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_Becker2007::SURFACE_TENSION_BOUNDARY, value);
+                }
+                else if (model->getSurfaceTensionMethod() == 2) // Akinci 2013
+                {
+                    SurfaceTension_Akinci2013* stMethod = static_cast<SurfaceTension_Akinci2013*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_Akinci2013::SURFACE_TENSION_BOUNDARY, value);
+                }
+                else if (model->getSurfaceTensionMethod() == 3) // He 2014
+                {
+                    SurfaceTension_He2014* stMethod = static_cast<SurfaceTension_He2014*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_He2014::SURFACE_TENSION_BOUNDARY, value);
+                }
+                else if (model->getSurfaceTensionMethod() == 4) // Jeske 2023
+                {
+                    SurfaceTension_Jeske2023* stMethod = static_cast<SurfaceTension_Jeske2023*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_Jeske2023::SURFACE_TENSION_BOUNDARY, value);
+                }
+                else {
+                    LOG_WARN << "Not a valid surface tension method.";
+                }
+            }
+            else if (paramName == "surfaceTensionViscosity")
+            {
+                if (model->getSurfaceTensionMethod() == 4) // Jeske 2023
+                {
+                    SurfaceTension_Jeske2023* stMethod = static_cast<SurfaceTension_Jeske2023*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_Jeske2023::VISCOSITY_COEFFICIENT, value);
+                }
+                else {
+                    LOG_WARN << "Not a valid surface tension method.";
+                }
+            }
+            else if (paramName == "surfaceTensionViscosityBoundary")
+            {
+                if (model->getSurfaceTensionMethod() == 4) // Jeske 2023
+                {
+                    SurfaceTension_Jeske2023* stMethod = static_cast<SurfaceTension_Jeske2023*>(model->getSurfaceTensionBase());
+                    stMethod->setValue(SurfaceTension_Jeske2023::VISCOSITY_COEFFICIENT_BOUNDARY, value);
+                }
+                else {
+                    LOG_WARN << "Not a valid surface tension method.";
+                }
+            }
+            else
+            {
+                LOG_WARN << "Unknown parameter: " << paramName; 
+            }
+        }
+    }
+}
diff --git a/SPlisHSPlasH/DynamicParameterSystem.h b/SPlisHSPlasH/DynamicParameterSystem.h
new file mode 100644
index 00000000..826a6f90
--- /dev/null
+++ b/SPlisHSPlasH/DynamicParameterSystem.h
@@ -0,0 +1,81 @@
+#ifndef __DynamicParameterSystem_h__
+#define __DynamicParameterSystem_h__
+
+#include "Common.h"
+#include "extern/tinyexpr/tinyexpr.h"
+#include <vector>
+#include <string>
+
+namespace SPH
+{
+    //class FluidModel;
+    //class ViscosityBase;
+
+    struct TinyExprDeleter {
+        void operator()(te_expr* expr) const {
+            if (expr) {
+                te_free(expr);
+            }
+        }
+    };
+
+    struct ParameterSchedule
+    {
+        std::string fluidId;
+        std::string parameter_name;
+        std::string expression;
+        std::vector<std::pair<Real, Real>> timeline;
+        int current_timeline_index;
+        bool use_expression;
+        bool use_step_function;
+        std::unique_ptr<te_expr, TinyExprDeleter> compiled_expr;
+        Real current_value;
+        Real default_value;
+        bool active;
+
+        ParameterSchedule() : compiled_expr(nullptr), current_value(0.0), default_value(0.0), use_expression(false), active(true), current_timeline_index(-1) {}
+        /*~ParameterSchedule() 
+        {
+            if (compiled_expr) 
+            {
+                te_free(compiled_expr);
+                compiled_expr = nullptr;
+            }
+        }*/
+    };
+
+    class DynamicParameterSystem
+    {
+    private:
+        std::vector<ParameterSchedule> m_schedules;
+
+        double m_t_double;
+        double m_dt_double;
+
+    public:
+        DynamicParameterSystem();
+        virtual ~DynamicParameterSystem();
+
+        void step();
+        void reset();
+
+        // Add schedules
+        bool addTimelineSchedule(const std::string& fluidId, const std::string& paramName, const std::vector<std::pair<Real, Real>>& timeline,
+            Real defaultValue = 1.0, bool stepFunction = false);
+
+        bool addExpressionSchedule(const std::string& fluidId, const std::string& paramName,
+            const std::string& expression, Real defaultValue = 1.0);
+
+        size_t numberOfSchedules() const { return m_schedules.size(); }
+        ParameterSchedule& getSchedule(const unsigned int i) { return m_schedules[i]; }
+        const ParameterSchedule& getSchedule(const unsigned int i) const { return m_schedules[i]; }
+
+    private:
+        bool compileExpression(ParameterSchedule& schedule);
+        Real interpolateTimeline(ParameterSchedule& schedule, Real currentTime);
+        Real stepTimeline(ParameterSchedule& schedule, Real currentTime);
+        void applyParameterValue(const std::string& fluidId, const std::string& paramName, Real value);
+    };
+}
+
+#endif
diff --git a/SPlisHSPlasH/ICSPH/TimeStepICSPH.cpp b/SPlisHSPlasH/ICSPH/TimeStepICSPH.cpp
index 3ff7e52f..9bba29c1 100644
--- a/SPlisHSPlasH/ICSPH/TimeStepICSPH.cpp
+++ b/SPlisHSPlasH/ICSPH/TimeStepICSPH.cpp
@@ -177,6 +177,7 @@ void TimeStepICSPH::step()
 
 	sim->emitParticles();
 	sim->animateParticles();
+	sim->updateDynamicParameters();
 
 	// Compute new time	
 	tm->setTime (tm->getTime () + h);
diff --git a/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp b/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp
index e3b10584..df151223 100644
--- a/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp
+++ b/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp
@@ -132,6 +132,7 @@ void TimeStepIISPH::step()
 
 	sim->emitParticles();
 	sim->animateParticles();
+	sim->updateDynamicParameters();
 
 	// Compute new time	
 	tm->setTime (tm->getTime () + h);
diff --git a/SPlisHSPlasH/PBF/TimeStepPBF.cpp b/SPlisHSPlasH/PBF/TimeStepPBF.cpp
index 0d4ccc92..65f265c2 100644
--- a/SPlisHSPlasH/PBF/TimeStepPBF.cpp
+++ b/SPlisHSPlasH/PBF/TimeStepPBF.cpp
@@ -181,6 +181,7 @@ void TimeStepPBF::step()
 
 	sim->emitParticles();
 	sim->animateParticles();
+	sim->updateDynamicParameters();
 
 	// Compute new time	
 	tm->setTime (tm->getTime () + h);
diff --git a/SPlisHSPlasH/PCISPH/TimeStepPCISPH.cpp b/SPlisHSPlasH/PCISPH/TimeStepPCISPH.cpp
index 11ce2bcc..3f0573aa 100644
--- a/SPlisHSPlasH/PCISPH/TimeStepPCISPH.cpp
+++ b/SPlisHSPlasH/PCISPH/TimeStepPCISPH.cpp
@@ -134,6 +134,7 @@ void TimeStepPCISPH::step()
 
 	sim->emitParticles();
 	sim->animateParticles();
+	sim->updateDynamicParameters();
 
 	// Compute new time		
 	tm->setTime(tm->getTime() + h);
diff --git a/SPlisHSPlasH/PF/TimeStepPF.cpp b/SPlisHSPlasH/PF/TimeStepPF.cpp
index 4bf2c94b..26f098ff 100644
--- a/SPlisHSPlasH/PF/TimeStepPF.cpp
+++ b/SPlisHSPlasH/PF/TimeStepPF.cpp
@@ -160,6 +160,7 @@ void TimeStepPF::step()
 	// update emitters
 	sim->emitParticles();
 	sim->animateParticles();
+	sim->updateDynamicParameters();
 	// Compute new time	
 	sim->updateTimeStepSize();
 	tm->setTime(tm->getTime () + tm->getTimeStepSize());
diff --git a/SPlisHSPlasH/Simulation.cpp b/SPlisHSPlasH/Simulation.cpp
index 8485a214..939738d1 100644
--- a/SPlisHSPlasH/Simulation.cpp
+++ b/SPlisHSPlasH/Simulation.cpp
@@ -85,6 +85,7 @@ Simulation::Simulation ()
 	m_counter = 0;
 
 	m_animationFieldSystem = new AnimationFieldSystem();
+	m_dynamicParameterSystem = nullptr;
 	m_boundaryHandlingMethod = static_cast<int>(BoundaryHandlingMethods::Bender2019);
 }
 
@@ -94,6 +95,7 @@ Simulation::~Simulation ()
 	delete m_debugTools;
 #endif
 	delete m_animationFieldSystem;
+	delete m_dynamicParameterSystem;
 	delete m_timeStep;
 	delete m_neighborhoodSearch;
 	delete TimeManager::getCurrent();
@@ -680,6 +682,12 @@ void Simulation::animateParticles()
 	STOP_TIMING_AVG
 }
 
+void Simulation::updateDynamicParameters()
+{
+	if (m_dynamicParameterSystem != nullptr)
+		m_dynamicParameterSystem->step();
+}
+
 void Simulation::addBoundaryModel(BoundaryModel *bm)
 {
 	m_boundaryModels.push_back(bm);
diff --git a/SPlisHSPlasH/Simulation.h b/SPlisHSPlasH/Simulation.h
index 1ee05db3..19bcb26e 100644
--- a/SPlisHSPlasH/Simulation.h
+++ b/SPlisHSPlasH/Simulation.h
@@ -8,6 +8,7 @@
 #include "NeighborhoodSearch.h"
 #include "BoundaryModel.h"
 #include "AnimationFieldSystem.h"
+#include "DynamicParameterSystem.h"
 #include "Utilities/FileSystem.h"
 #ifdef USE_DEBUG_TOOLS
 #include "SPlisHSPlasH/Utilities/DebugTools.h"
@@ -289,6 +290,7 @@ namespace SPH
 		std::vector<FluidInfo> m_fluidInfos;
 		NeighborhoodSearch *m_neighborhoodSearch;
 		AnimationFieldSystem *m_animationFieldSystem;
+		DynamicParameterSystem *m_dynamicParameterSystem;
 		int m_cflMethod;
 		Real m_cflFactor;
 		Real m_cflMinTimeStepSize;
@@ -365,6 +367,10 @@ namespace SPH
 		FluidInfo& getFluidInfo(const unsigned int i) { return m_fluidInfos[i]; }
 
 		AnimationFieldSystem* getAnimationFieldSystem() { return m_animationFieldSystem; }
+
+		DynamicParameterSystem* getDynamicParameterSystem() { return m_dynamicParameterSystem; }
+		void setDynamicParameterSystem(DynamicParameterSystem* dps) { m_dynamicParameterSystem = dps; }
+		void updateDynamicParameters();
 		
 		BoundaryHandlingMethods getBoundaryHandlingMethod() const { return (BoundaryHandlingMethods) m_boundaryHandlingMethod; }
 		void setBoundaryHandlingMethod(BoundaryHandlingMethods val) { m_boundaryHandlingMethod = (int) val; }
diff --git a/SPlisHSPlasH/Utilities/SceneLoader.cpp b/SPlisHSPlasH/Utilities/SceneLoader.cpp
index cc2f5bca..70f6064c 100644
--- a/SPlisHSPlasH/Utilities/SceneLoader.cpp
+++ b/SPlisHSPlasH/Utilities/SceneLoader.cpp
@@ -133,6 +133,21 @@ void SceneLoader::readScene(const char *fileName, Scene &scene)
 		}
 	}
 
+	//////////////////////////////////////////////////////////////////////////
+	// read dynamic parameters
+	//////////////////////////////////////////////////////////////////////////
+	if (m_jsonData.find("DynamicParameters") != m_jsonData.end())
+	{
+		nlohmann::json dynParams = m_jsonData["DynamicParameters"];
+		for (auto& param : dynParams)
+		{
+			DynamicParameterObject* data = new DynamicParameterObject();
+			data->initParameters();
+			readParameterObject(param, data);
+			scene.dynamicParameters.push_back(data);
+		}
+	}
+
 	//////////////////////////////////////////////////////////////////////////
 	// read materials
 	//////////////////////////////////////////////////////////////////////////
diff --git a/SPlisHSPlasH/Utilities/SceneLoader.h b/SPlisHSPlasH/Utilities/SceneLoader.h
index ebc445c6..58891984 100644
--- a/SPlisHSPlasH/Utilities/SceneLoader.h
+++ b/SPlisHSPlasH/Utilities/SceneLoader.h
@@ -29,6 +29,7 @@ namespace Utilities
 			std::vector<FluidBlockParameterObject*> fluidBlocks;
 			std::vector<EmitterParameterObject*> emitters;
 			std::vector<AnimationFieldParameterObject*> animatedFields;
+			std::vector<DynamicParameterObject*> dynamicParameters;
 			std::vector<MaterialParameterObject*> materials;
 			Real particleRadius;
 			bool sim2D;
diff --git a/SPlisHSPlasH/Utilities/SceneParameterObjects.cpp b/SPlisHSPlasH/Utilities/SceneParameterObjects.cpp
index 6be6bbc2..5839661c 100644
--- a/SPlisHSPlasH/Utilities/SceneParameterObjects.cpp
+++ b/SPlisHSPlasH/Utilities/SceneParameterObjects.cpp
@@ -246,6 +246,48 @@ void AnimationFieldParameterObject::initParameters()
 	setDescription(ANIMATIONFIELD_ENDTIME, "End time of the animation field.");	
 }
 
+//////////////////////////////////////////////////////////////////////////
+// DynamicParameterObject
+//////////////////////////////////////////////////////////////////////////
+int DynamicParameterObject::DYNPARAM_FLUID = -1;
+int DynamicParameterObject::DYNPARAM_NAME = -1;
+int DynamicParameterObject::DYNPARAM_EXPRESSION = -1;
+int DynamicParameterObject::DYNPARAM_DEFAULT_VALUE = -1;
+int DynamicParameterObject::DYNPARAM_STEP_FUNCTION = -1;
+int DynamicParameterObject::DYNPARAM_TIMELINE_TIMES = -1;
+int DynamicParameterObject::DYNPARAM_TIMELINE_VALUES = -1;
+
+void DynamicParameterObject::initParameters()
+{
+	DYNPARAM_FLUID = createStringParameter("fluidId", "Fluid ID", &fluidId);
+	setGroup(DYNPARAM_FLUID, "Dynamic Parameter");
+	setDescription(DYNPARAM_FLUID, "Fluid on which to apply the parameter change.");
+
+	DYNPARAM_NAME = createStringParameter("parameterName", "Parameter name", &parameterName);
+	setGroup(DYNPARAM_NAME, "Dynamic Parameter");
+	setDescription(DYNPARAM_NAME, "Name of the parameter to control (e.g., stiffness, viscosity, density0).");
+
+	DYNPARAM_EXPRESSION = createStringParameter("expression", "Expression", &expression);
+	setGroup(DYNPARAM_EXPRESSION, "Dynamic Parameter");
+	setDescription(DYNPARAM_EXPRESSION, "Mathematical expression for parameter control (using t for time, dt for timestep).");
+
+	DYNPARAM_DEFAULT_VALUE = createNumericParameter<Real>("defaultValue", "Default value", &defaultValue);
+	setGroup(DYNPARAM_DEFAULT_VALUE, "Dynamic Parameter");
+	setDescription(DYNPARAM_DEFAULT_VALUE, "Default value of the parameter.");
+
+	DYNPARAM_STEP_FUNCTION = createBoolParameter("stepFunction", "Step function", &stepFunction);
+	setGroup(DYNPARAM_STEP_FUNCTION, "Dynamic Parameter");
+	setDescription(DYNPARAM_STEP_FUNCTION, "Use step function instead of linear interpolation for timeline.");
+
+	DYNPARAM_TIMELINE_TIMES = createStringParameter("timelineTimes", "Timeline times", &timelineTimes);
+	setGroup(DYNPARAM_TIMELINE_TIMES, "Dynamic Parameter");
+	setDescription(DYNPARAM_TIMELINE_TIMES, "Time points for timeline-based parameter changes.");
+
+	DYNPARAM_TIMELINE_VALUES = createStringParameter("timelineValues", "Timeline values", &timelineValues);
+	setGroup(DYNPARAM_TIMELINE_VALUES, "Dynamic Parameter");
+	setDescription(DYNPARAM_TIMELINE_VALUES, "Values corresponding to timeline times.");
+}
+
 //////////////////////////////////////////////////////////////////////////
 // MaterialParameterObject
 //////////////////////////////////////////////////////////////////////////
diff --git a/SPlisHSPlasH/Utilities/SceneParameterObjects.h b/SPlisHSPlasH/Utilities/SceneParameterObjects.h
index 1eea317a..38e23112 100644
--- a/SPlisHSPlasH/Utilities/SceneParameterObjects.h
+++ b/SPlisHSPlasH/Utilities/SceneParameterObjects.h
@@ -260,6 +260,54 @@ namespace Utilities
 		virtual void initParameters();
 	};
 
+	/** \brief class to store a dynamic parameter object
+	 */
+	class DynamicParameterObject : public GenParam::ParameterObject
+	{
+	public:
+		std::string fluidId;
+		std::string parameterName;
+		std::string expression;
+		Real defaultValue;
+		bool stepFunction;
+		std::string timelineTimes;
+		std::string timelineValues;
+
+		DynamicParameterObject()
+		{
+			// Default values
+			fluidId = "";
+			parameterName = "";
+			expression = "";
+			defaultValue = 1.0;
+			stepFunction = false;
+			timelineTimes = "";
+			timelineValues = "";
+		}
+
+		DynamicParameterObject(std::string fluidId_, std::string parameterName_, std::string expression_, Real defaultValue_, bool stepFunction_,
+			std::string timelineTimes_, std::string timelineValues_)
+		{
+			fluidId = fluidId_;
+			parameterName = parameterName_;
+			expression = expression_;
+			defaultValue = defaultValue_;
+			stepFunction = stepFunction_;
+			timelineTimes = timelineTimes_;
+			timelineValues = timelineValues_;
+		}
+
+		static int DYNPARAM_FLUID;
+		static int DYNPARAM_NAME;
+		static int DYNPARAM_EXPRESSION;
+		static int DYNPARAM_DEFAULT_VALUE;
+		static int DYNPARAM_STEP_FUNCTION;
+		static int DYNPARAM_TIMELINE_TIMES;
+		static int DYNPARAM_TIMELINE_VALUES;
+
+		virtual void initParameters();
+	};
+
 	/** \brief Class to store particle coloring information */
 	class MaterialParameterObject : public GenParam::ParameterObject
 	{
diff --git a/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp b/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp
index ab3bcaca..dad5c3a4 100644
--- a/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp
+++ b/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp
@@ -135,6 +135,7 @@ void TimeStepWCSPH::step()
 
 	sim->emitParticles();
 	sim->animateParticles();
+	sim->updateDynamicParameters();
 
 	// Compute new time	
 	tm->setTime (tm->getTime () + h);
diff --git a/Simulator/SceneConfiguration.cpp b/Simulator/SceneConfiguration.cpp
index 191731c4..e6895034 100644
--- a/Simulator/SceneConfiguration.cpp
+++ b/Simulator/SceneConfiguration.cpp
@@ -36,6 +36,10 @@ SceneConfiguration::~SceneConfiguration ()
 		delete m_scene.animatedFields[i];
 	m_scene.animatedFields.clear();
 
+	for (unsigned int i = 0; i < m_scene.dynamicParameters.size(); i++)
+		delete m_scene.dynamicParameters[i];
+	m_scene.dynamicParameters.clear();
+
 	m_current = nullptr;
 }
 
diff --git a/Simulator/SimulatorBase.cpp b/Simulator/SimulatorBase.cpp
index 5f1317c7..773f721a 100644
--- a/Simulator/SimulatorBase.cpp
+++ b/Simulator/SimulatorBase.cpp
@@ -7,6 +7,7 @@
 #include "Utilities/PartioReaderWriter.h"
 #include "SPlisHSPlasH/Emitter.h"
 #include "SPlisHSPlasH/EmitterSystem.h"
+#include "SPlisHSPlasH/DynamicParameterSystem.h"
 #include "SPlisHSPlasH/Simulation.h"
 #include "SPlisHSPlasH/NonPressureForceBase.h"
 #include "NumericParameter.h"
@@ -811,6 +812,7 @@ void SimulatorBase::buildModel()
 
 	createEmitters();
 	createAnimationFields();
+	createDynamicParameters();
 
 	Simulation *sim = Simulation::getCurrent();
 
@@ -1414,6 +1416,89 @@ void SimulatorBase::createAnimationFields()
 	}
 }
 
+void SimulatorBase::createDynamicParameters()
+{
+	Simulation* sim = Simulation::getCurrent();
+	const Utilities::SceneLoader::Scene& scene = SceneConfiguration::getCurrent()->getScene();
+
+	//////////////////////////////////////////////////////////////////////////
+	// Dynamic parameters
+	//////////////////////////////////////////////////////////////////////////
+
+	if (scene.dynamicParameters.empty())
+		return;
+
+	if (sim->getDynamicParameterSystem() == nullptr)
+	{
+		DynamicParameterSystem* dps = new DynamicParameterSystem();
+		sim->setDynamicParameterSystem(dps);
+	}
+
+	DynamicParameterSystem* dps = sim->getDynamicParameterSystem();
+	for (unsigned int i = 0; i < scene.dynamicParameters.size(); i++)
+	{
+		DynamicParameterObject* data = scene.dynamicParameters[i];
+
+		const std::string& fluidId = data->fluidId;
+		const std::string& paramName = data->parameterName;
+		Real defaultValue = data->defaultValue;
+
+		if (!data->expression.empty())
+		{
+			dps->addExpressionSchedule(fluidId, paramName, data->expression, defaultValue);
+		}
+		else if (!data->timelineTimes.empty() && !data->timelineValues.empty())
+		{
+			std::vector<Real> times = parseCommaSeparatedString(data->timelineTimes);
+			std::vector<Real> values = parseCommaSeparatedString(data->timelineValues);
+
+			if (times.size() != values.size()) {
+				LOG_WARN << "Timeline times and values have different sizes for parameter: " << paramName;
+				continue;
+			}
+
+			std::vector<std::pair<Real, Real>> timeline;
+			for (size_t i = 0; i < times.size(); ++i) {
+				timeline.push_back({ times[i], values[i] });
+			}
+
+			dps->addTimelineSchedule(fluidId, paramName, timeline, defaultValue, data->stepFunction);
+		}
+	}
+
+}
+
+std::vector<Real> SimulatorBase::parseCommaSeparatedString(const std::string& str)
+{
+	std::vector<Real> result;
+	if (str.empty()) return result;
+
+	std::stringstream ss(str);
+	std::string item;
+
+	while (std::getline(ss, item, ',')) {
+		// Remove spaces
+		item.erase(item.begin(), std::find_if(item.begin(), item.end(), [](unsigned char ch) {
+			return !std::isspace(ch);
+			}));
+		item.erase(std::find_if(item.rbegin(), item.rend(), [](unsigned char ch) {
+			return !std::isspace(ch);
+			}).base(), item.end());
+
+		if (!item.empty()) {
+			try {
+				Real value = std::stod(item);
+				result.push_back(value);
+			}
+			catch (const std::exception& e) {
+				LOG_WARN << "Could not parse value: " << item;
+			}
+		}
+	}
+
+	return result;
+}
+
 
 void SimulatorBase::createFluidBlocks(std::map<std::string, unsigned int> &fluidIDs, std::vector<std::vector<Vector3r>> &fluidParticles, std::vector<std::vector<Vector3r>> &fluidVelocities, std::vector<std::vector<unsigned int>> &fluidObjectIds)
 {
diff --git a/Simulator/SimulatorBase.h b/Simulator/SimulatorBase.h
index a1ed2031..479cbe5d 100644
--- a/Simulator/SimulatorBase.h
+++ b/Simulator/SimulatorBase.h
@@ -101,6 +101,8 @@ namespace SPH
 		void createFluidBlocks(std::map<std::string, unsigned int> &fluidIDs, std::vector<std::vector<Vector3r>> &fluidParticles, std::vector<std::vector<Vector3r>> &fluidVelocities, std::vector<std::vector<unsigned int>> &fluidObjectIds);
 		void createEmitters();
 		void createAnimationFields();
+		void createDynamicParameters();
+		std::vector<Real> parseCommaSeparatedString(const std::string& str); // To parse the dynamicParameters timelines
 		void buildModel();
 		void setCommandLineParameter();
 		void setCommandLineParameter(GenParam::ParameterObject *paramObj);