diff --git a/src/discretedomain.jl b/src/discretedomain.jl
index c7f90007c5..67bef42c9e 100644
--- a/src/discretedomain.jl
+++ b/src/discretedomain.jl
@@ -85,7 +85,7 @@ $(TYPEDEF)
 Represents a sample operator. A discrete-time signal is created by sampling a continuous-time signal.
 
 # Constructors
-`Sample(clock::Union{TimeDomain, InferredTimeDomain} = InferredDiscrete)`
+`Sample(clock::Union{TimeDomain, InferredTimeDomain} = InferredDiscrete())`
 `Sample(dt::Real)`
 
 `Sample(x::Num)`, with a single argument, is shorthand for `Sample()(x)`.
@@ -106,7 +106,7 @@ julia> Δ = Sample(0.01)
 """
 struct Sample <: Operator
     clock::Any
-    Sample(clock::Union{TimeDomain, InferredTimeDomain} = InferredDiscrete) = new(clock)
+    Sample(clock::Union{TimeDomain, InferredTimeDomain} = InferredDiscrete()) = new(clock)
 end
 
 function Sample(arg::Real)
@@ -190,7 +190,7 @@ struct ShiftIndex
     clock::Union{InferredTimeDomain, TimeDomain, IntegerSequence}
     steps::Int
     function ShiftIndex(
-            clock::Union{TimeDomain, InferredTimeDomain, IntegerSequence} = Inferred, steps::Int = 0)
+            clock::Union{TimeDomain, InferredTimeDomain, IntegerSequence} = Inferred(), steps::Int = 0)
         new(clock, steps)
     end
     ShiftIndex(dt::Real, steps::Int = 0) = new(Clock(dt), steps)
@@ -232,7 +232,7 @@ function input_timedomain(s::Shift, arg = nothing)
     if has_time_domain(arg)
         return get_time_domain(arg)
     end
-    InferredDiscrete
+    InferredDiscrete()
 end
 
 """
@@ -244,7 +244,7 @@ function output_timedomain(s::Shift, arg = nothing)
     if has_time_domain(t, arg)
         return get_time_domain(t, arg)
     end
-    InferredDiscrete
+    InferredDiscrete()
 end
 
 input_timedomain(::Sample, _ = nothing) = Continuous()
@@ -254,7 +254,7 @@ function input_timedomain(h::Hold, arg = nothing)
     if has_time_domain(arg)
         return get_time_domain(arg)
     end
-    InferredDiscrete # the Hold accepts any discrete
+    InferredDiscrete() # the Hold accepts any discrete
 end
 output_timedomain(::Hold, _ = nothing) = Continuous()
 
diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl
index ad3fbfdeef..8b893b6cd2 100644
--- a/src/systems/nonlinear/nonlinearsystem.jl
+++ b/src/systems/nonlinear/nonlinearsystem.jl
@@ -219,6 +219,12 @@ function NonlinearSystem(eqs; kwargs...)
                 push!(new_ps, p)
             end
         else
+            if symbolic_type(p) == ArraySymbolic() &&
+               Symbolics.shape(unwrap(p)) != Symbolics.Unknown()
+                for i in eachindex(p)
+                    delete!(new_ps, p[i])
+                end
+            end
             push!(new_ps, p)
         end
     end
diff --git a/test/nonlinearsystem.jl b/test/nonlinearsystem.jl
index c6318b1314..c8dda530e2 100644
--- a/test/nonlinearsystem.jl
+++ b/test/nonlinearsystem.jl
@@ -4,6 +4,7 @@ using DiffEqBase, SparseArrays
 using Test
 using NonlinearSolve
 using ForwardDiff
+using SymbolicIndexingInterface
 using ModelingToolkit: value
 using ModelingToolkit: get_default_or_guess, MTKParameters
 
@@ -380,3 +381,12 @@ end
     @test_throws ["single equation", "unknown"] IntervalNonlinearFunctionExpr(
         sys, (0.0, 1.0))
 end
+
+@testset "Vector parameter used unscalarized and partially scalarized" begin
+    @variables x y
+    @parameters p[1:2] (f::Function)(..)
+
+    @mtkbuild sys = NonlinearSystem([x^2 - p[1]^2 ~ 0, y^2 ~ f(p)])
+    @test !any(isequal(p[1]), parameters(sys))
+    @test is_parameter(sys, p)
+end
diff --git a/test/scc_nonlinear_problem.jl b/test/scc_nonlinear_problem.jl
index 92030a84aa..b2b326d090 100644
--- a/test/scc_nonlinear_problem.jl
+++ b/test/scc_nonlinear_problem.jl
@@ -282,3 +282,12 @@ end
     sccprob = SCCNonlinearProblem(fullsys, u0, p)
     @test isequal(parameters(fullsys), parameters(sccprob.f.sys))
 end
+
+@testset "Vector parameters in function arguments" begin
+    @variables x y
+    @parameters p[1:2] (f::Function)(..)
+
+    @mtkbuild sys = NonlinearSystem([x^2 - p[1]^2 ~ 0, y^2 ~ f(p)])
+    prob = SCCNonlinearProblem(sys, [x => 1.0, y => 1.0], [p => ones(2), f => sum])
+    @test_nowarn solve(prob, NewtonRaphson())
+end