Questions about calculating uncertainty and correlation matrix of model parameters from optimization

I am running a nonlinear earth system model to optimize 42 parameters p with 7 different kinds of observations $O_j$ where j=1,...,7. Each observations (or modelled estimates) have different variances and mean values, and have different time series length. The longest observation is a time series from 2005 to 2010 (say total there is 2000 data points). The shortest only have one data point. I am giving these 7 observation streams equal weights. I am wondering how can I calculate the Jacobians and corresponding parameter uncertainty&correlation matrix?

I checked there is a similar question: Working out error on fit parameters for nonlinear fit. But I am wondering is there any strict proof/derivation of this method? Besides, what is the dimension of Jacobian J?

Suppose I am using cost function for the 6 observations like (written in Julia)

function metric(y::AbstractArray, yσ::AbstractArray, ŷ::AbstractArray, ::NSE)
    NSE = one(eltype(ŷ)) .- sum(abs2.((y .- ŷ))) / sum(abs2.((y .- mean(y))))
    return NSE
end

function metric(y::AbstractArray, yσ::AbstractArray, ŷ::AbstractArray, ::NSEInv)
    NSEInv = one(eltype(ŷ)) - metric(y, yσ, ŷ, NSE())
    return NSEInv
end

and a different cost function for the shortest observation time series:

NAME1R = \\frac{(|μ_ŷ - μ_y|)}{1 + μ_y}

I have following questions:

1. Then J should be of dimension of [num_param, num_total_time_point]? By num_total_time_point I mean all time series data points of all 7 observations. and should be $\frac{dŷ}{dp_i}$ for $i-th$ parameter? or it is $\frac{dŷ}{dp_i}*\frac{1}{sum(abs2.((y .- mean(y)))}$ for those 6 observations and $\frac{dŷ}{dp_i}*\frac{1}{1 + μ_y}$ for the one shortest one?

2. J is the jacobian of residuals to parameters? or model outputs? i.e. it is $\frac{dŷ}{dp_i}$ or $\frac{d(ŷ-y)}{dp_i}$? at a specific time point? or $\frac{d(ŷ-y)^2}{dp_i}$ ?

3. After getting J, then I can use J' * J to approximate Covariance Matrix? not sure whether I can directly calculate Hessian matrix? because ForwardDiff.jl from Julia has this function...

Not sure about the choice about 1...I found a similar paper (see equation A6 in the appendix of https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-2486.2001.00434.x), however, this one does not consider the variance of observations in the matrix $Q$...

Any suggestions would be appreciated!!!