diff --git a/.github/sync.yml b/.github/sync.yml
deleted file mode 100644
index e8d5dc0..0000000
--- a/.github/sync.yml
+++ /dev/null
@@ -1,12 +0,0 @@
-group:
-  repos: |
-    ODINN-SciML/Huginn.jl
-    ODINN-SciML/Muninn.jl
-    ODINN-SciML/ODINN.jl
-  files: 
-    - source: environment.yml
-      dest: environment.yml
-    - source: .gitignore
-      dest: .gitignore
-    - source: .github/workflows/CI.yml
-      dest: .github/workflows/CI.yml
diff --git a/Project.toml b/Project.toml
index dd7b346..c24376f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
+DateFormats = "44557152-fe0a-4de1-8405-416d90313ce6"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
@@ -17,6 +18,7 @@ Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
+NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
diff --git a/src/Sleipnir.jl b/src/Sleipnir.jl
index 5681457..2e56342 100644
--- a/src/Sleipnir.jl
+++ b/src/Sleipnir.jl
@@ -23,6 +23,8 @@ using Tar
 import NCDatasets
 using Unitful: m, rad, °
 using CoordRefSystems
+using Dates, DateFormats
+using NetCDF
 
 # ##############################################
 # ############    PARAMETERS     ###############
@@ -39,11 +41,13 @@ const global prepro_dir::String = joinpath(homedir(), ".ODINN", "ODINN_prepro")
 include("setup/config.jl")
 # All parameters needed for the models
 include("parameters/Parameters.jl")
-# Anything related to managing glacier topographical and climate data
+# Anything related to managing glacier topographical and climate variables
 include("glaciers/glacier/Glacier.jl")
+# Anything related to managing glacier data used for data assimilation
+include("glaciers/data/Data.jl")
 # All structures and functions related to ODINN models
 include("models/Model.jl")
 # Everything related to running simulations in ODINN
 include("simulations/Simulation.jl")
 
-end # module
+end # module
\ No newline at end of file
diff --git a/src/glaciers/data/Data.jl b/src/glaciers/data/Data.jl
new file mode 100644
index 0000000..e2d1ae7
--- /dev/null
+++ b/src/glaciers/data/Data.jl
@@ -0,0 +1,5 @@
+export AbstractData
+
+abstract type AbstractData end
+
+include("SurfaceVelocityData.jl")
\ No newline at end of file
diff --git a/src/glaciers/data/SurfaceVelocityData.jl b/src/glaciers/data/SurfaceVelocityData.jl
new file mode 100644
index 0000000..61987d5
--- /dev/null
+++ b/src/glaciers/data/SurfaceVelocityData.jl
@@ -0,0 +1,75 @@
+export SurfaceVelocityData
+
+mutable struct SurfaceVelocityData{F <: AbstractFloat} <: AbstractData
+    x::Union{Vector{F}, Nothing}
+    y::Union{Vector{F}, Nothing}
+    lat::Union{Vector{F}, Nothing}
+    lon::Union{Vector{F}, Nothing}
+    vx::Union{Array{F, 3}, Nothing}
+    vy::Union{Array{F, 3}, Nothing}
+    vabs::Union{Array{F, 3}, Nothing}
+    vx_error::Union{Array{F, 1}, Nothing}
+    vy_error::Union{Array{F, 1}, Nothing}
+    vabs_error::Union{Array{F, 1}, Nothing}
+    date::Union{Vector{DateTime}, Nothing}
+    date1::Union{Vector{DateTime}, Nothing}
+    date2::Union{Vector{DateTime}, Nothing}
+    date_error::Union{Vector{Day}, Vector{Millisecond}, Nothing}
+end
+
+"""
+function SurfaceVelocityData(;
+    x::Union{Vector{F}, Nothing} = nothing,
+    y::Union{Vector{F}, Nothing} = nothing, 
+    lat::Union{Vector{F}, Nothing} = nothing,
+    lon::Union{Vector{F}, Nothing} = nothing, 
+    vx::Union{Array{F, 3}, Nothing} = nothing,
+    vy::Union{Array{F, 3}, Nothing} = nothing,
+    vabs::Union{Array{F, 3}, Nothing} = nothing,
+    vx_error::Union{Array{F, 1}, Nothing} = nothing,
+    vy_error::Union{Array{F, 1}, Nothing} = nothing,
+    vabs_error::Union{Array{F, 1}, Nothing} = nothing,
+    date::Union{Vector{DateTime}, Nothing} = nothing,
+    date1::Union{Vector{DateTime}, Nothing} = nothing,
+    date2::Union{Vector{DateTime}, Nothing} = nothing,
+    date_error::Union{Vector{Day}, Vector{Millisecond}, Nothing} = nothing,
+    ) where {F <: AbstractFloat}
+
+Constructor for ice surface velocity data based on Rabatel et. al (2023). 
+
+
+Important remarks:
+- Projections in longitude and latitude assume we are working in the north hemisphere. 
+  If working with south hemisphere glaciers, this needs to be changed.
+- The error in velocity is unique per timestamp, rather than being pixel distributed. 
+- The error in the absolute velocities `vabs_error` is overestimated.
+
+References:
+    - Rabatel, A., Ducasse, E., Millan, R. & Mouginot, J. 
+    Satellite-Derived Annual Glacier Surface Flow Velocity Products for the European Alps, 2015–2021. 
+    Data 8, 66 (2023).
+"""
+function SurfaceVelocityData(;
+    x::Union{Vector{F}, Nothing} = nothing,
+    y::Union{Vector{F}, Nothing} = nothing, 
+    lat::Union{Vector{F}, Nothing} = nothing,
+    lon::Union{Vector{F}, Nothing} = nothing, 
+    vx::Union{Array{F, 3}, Nothing} = nothing,
+    vy::Union{Array{F, 3}, Nothing} = nothing,
+    vabs::Union{Array{F, 3}, Nothing} = nothing,
+    vx_error::Union{Array{F, 1}, Nothing} = nothing,
+    vy_error::Union{Array{F, 1}, Nothing} = nothing,
+    vabs_error::Union{Array{F, 1}, Nothing} = nothing,
+    date::Union{Vector{DateTime}, Nothing} = nothing,
+    date1::Union{Vector{DateTime}, Nothing} = nothing,
+    date2::Union{Vector{DateTime}, Nothing} = nothing,
+    date_error::Union{Vector{Day}, Vector{Millisecond}, Nothing} = nothing,
+    ) where {F <: AbstractFloat}
+
+    ft = typeof(x[begin])
+
+    return SurfaceVelocityData{ft}(x, y, lat, lon, vx, vy, vabs, vx_error, vy_error, vabs_error, date, date1, date2, date_error)
+end
+
+
+include("surfacevelocitydata_utils.jl")
\ No newline at end of file
diff --git a/src/glaciers/data/surfacevelocitydata_utils.jl b/src/glaciers/data/surfacevelocitydata_utils.jl
new file mode 100644
index 0000000..e1b721e
--- /dev/null
+++ b/src/glaciers/data/surfacevelocitydata_utils.jl
@@ -0,0 +1,109 @@
+export initialize_surfacevelocitydata
+export plot_timeseries, plot_count
+
+"""
+    initialize_surfacevelocitydata(file::String; interp=false)
+
+Initialize SurfaceVelocityData from Rabatel et. al (2023). 
+
+Arguments
+=================
+    - `file`: name of netcdf file with data
+    - `interp`: boolean variable indicating if we use the inporpolated data or not
+"""
+function initialize_surfacevelocitydata(file::String; interp=false, compute_vabs_error::Bool=true)
+
+    # Date of first adquisition
+    date1 = ncread(file, "date1")
+    # Date of second adquisition
+    date2 = ncread(file, "date2")
+
+    # Middle date
+    if !interp 
+        date_mean = 0.5 .* date1 .+ 0.5 .* date2
+        date_mean_offset = date1_offset = date2_offset = 0
+    else 
+        date_mean = ncread(file, "mid_date")
+        # Extract string from metadata of when days start counting
+        # Default format of attribute "units" is "days since 2015-08-14 00:00:00"
+        date_mean_since = ncgetatt(file, "mid_date", "units")[12:21] # e.g., "2015-08-14"
+        date1_offset_since = ncgetatt(file, "date1", "units")[12:21] # e.g., "2015-07-30"
+        date2_offset_since = ncgetatt(file, "date2", "units")[12:21] # e.g., "2015-08-29"
+        # Convertion to Julia datetime
+        date_mean_offset = datetime2julian(DateTime(date_mean_since)) - 2400000.5
+        date1_offset = datetime2julian(DateTime(date1_offset_since)) - 2400000.5
+        date2_offset = datetime2julian(DateTime(date_mean_offset)) - 2400000.5
+    end
+
+    # Convert dates from Modified Julian Days to datetipes
+    date1 = mjd.(date1 .+ date1_offset)
+    date2 = mjd.(date2 .+ date2_offset)
+    date_mean = mjd.(date_mean .+ date_mean_offset)
+    date_error = date2 .- date1
+
+    # We convert from Date to DateTime
+    date_mean = DateTime.(date_mean)
+    date1 = DateTime.(date1)
+    date2 = DateTime.(date2)
+
+    # Read data from netcdf file
+    x = ncread(file, "x")
+    y = ncread(file, "y")
+
+    # Velocity in the x direction (m/yr)
+    vx = ncread(file, "vx")              
+    vy = ncread(file, "vy")    
+
+    # Compute absolute velocity
+    vabs = (vx.^2 .+ vy.^2).^0.5  
+    # The sqrt operation in Julia promotes Float32 to Float64. We convert manually 
+    # to keep type consistency     
+    vabs = convert(typeof(vx), vabs)   
+    
+    # Error is reported once per timespan, so upper bounds are given by absolute error
+    if !interp
+        vx_error = ncread(file, "error_vx")
+        vy_error = ncread(file, "error_vy")
+        # Absolute error uncertanty using propagation of uncertanties 
+        if compute_vabs_error
+            vx_ratio_max = map(i -> max_or_empty(abs.(vx[:,:,i][vabs[:,:,i] .> 0.0]) ./ vabs[:,:,i][vabs[:,:,i] .> 0.0]), 1:size(vx)[3])
+            vy_ratio_max = map(i -> max_or_empty(abs.(vy[:,:,i][vabs[:,:,i] .> 0.0]) ./ vabs[:,:,i][vabs[:,:,i] .> 0.0]), 1:size(vy)[3])
+            vabs_error = ((vx_ratio_max .* vx_error).^2 .+ (vy_ratio_max .* vy_error).^2).^0.5
+            vabs_error = convert(typeof(vx_error), vabs_error)  
+        else
+            vabs_error = nothing
+        end
+    else         
+        vx_error = nothing
+        vy_error = nothing
+        vabs_error = nothing
+    end
+
+    # Run some basic tests 
+    nx, ny, ntimes = size(vx)
+    @assert length(date1) == length(date2) == ntimes
+    @assert nx == ny == 250
+
+    # Spatial preprocessing
+    proj_zone = ncgetatt(file, "mapping", "utm_zone_number")
+    transform(X,Y) = Sleipnir.UTMercator(X, Y; zone=proj_zone, hemisphere=:north) 
+    Coordinates = transform.(x, y)
+    latitudes = map(x -> x.lat.val, Coordinates)
+    longitudes = map(x -> x.lon.val, Coordinates)
+
+    return SurfaceVelocityData(x=x, y=y, lat=latitudes, lon=longitudes, vx=vx, vy=vy, vabs=vabs, vx_error=vx_error, vy_error=vy_error, vabs_error=vabs_error, date=date_mean, date1=date1, date2=date2, date_error=date_error)
+end
+
+"""
+    max_or_empyt(A::Array)
+
+Return maximum value for non-empty arrays. 
+This is just required to compute the error in the absolute velocity.
+"""
+function max_or_empty(A::Array)
+    if length(A) == 0
+        return 0.0
+    else
+        return maximum(A)
+    end
+end
\ No newline at end of file
diff --git a/src/parameters/PhysicalParameters.jl b/src/parameters/PhysicalParameters.jl
index f7c6784..796e975 100644
--- a/src/parameters/PhysicalParameters.jl
+++ b/src/parameters/PhysicalParameters.jl
@@ -13,28 +13,29 @@ end
 
 """
     PhysicalParameters(;
-        ρ::Float64 = 900.0,
-        g::Float64 = 9.81,
-        ϵ::Float64 = 1e-3,
-        η₀::F = 1.0, 
-        maxA::Float64 = 8e-17,
-        minA::Float64 = 8.5e-20,
-        maxTlaw::Float64 = 1.0,
-        minTlaw::Float64 = -25.0,
-        noise_A_magnitude::Float64 = 5e-18
-        )
+            ρ::F = 900.0,
+            g::F = 9.81,
+            ϵ::F = 1e-3,
+            η₀::F = 1.0, 
+            maxA::F = 8e-17,
+            minA::F = 8.5e-20,
+            maxTlaw::F = 1.0,
+            minTlaw::F = -25.0,
+            noise_A_magnitude::F = 5e-18
+            )
+
 Initialize the physical parameters of a model.
 Keyword arguments
 =================
     - `ρ`: Ice density
     - `g`: Gravitational constant
-    - `n`: Glen's exponent
-    - `A`: Glen's coefficient
     - `ϵ`: Small number
-    - `C`: Sliding coefficient
     - `η₀`:  
     - `maxA`: Maximum value for `A` (Glen's coefficient)
     - `minA`: Minimum value for `A` (Glen's coefficient)
+    - `maxTlaw`: Maximum value of Temperature used in simulations on fake law
+    - `minTlaw`: Minimum value of Temperature used in simulations on fake law
+    - `noise_A_magnitude`: Magnitude of noise added to A
 """
 function PhysicalParameters(;
             ρ::F = 900.0,