1. Calibration Site Parameters

Calibrate PMLV2 model using FLUXNET data

US Twt (CRO, 2009)

using PenmanMonteithLeuning, Ipaper

df_out, df, _par = deserialize(file_FLUXNET_CRO_USTwt)
(df_out = 275×10 DataFrame
 Row │ GPP_sim   ET_sim    Es         Ecr         Eca        Ei          Es_eq ⋯
 Float64   Float64   Float64    Float64     Float64    Float64     Float ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ 0.797755  0.309498  0.234323   0.00472006  0.0658945  0.00456045  0.234 ⋯
   2 │ 1.49914   0.389925  0.172465   0.0117285   0.205732   0.0         0.265
   3 │ 1.66223   0.604574  0.372541   0.0411049   0.17275    0.0181786   0.372
   4 │ 2.54749   0.795978  0.430526   0.0593096   0.297933   0.00820882  0.430
   5 │ 2.87331   0.881985  0.441199   0.0947368   0.300108   0.0459409   0.441 ⋯
   6 │ 2.67144   0.804064  0.355279   0.0605682   0.175928   0.212289    0.355
   7 │ 3.77559   1.03439   0.457163   0.177544    0.296867   0.102812    0.457
   8 │ 4.18786   1.07175   0.421187   0.204591    0.290327   0.155647    0.421
   9 │ 6.01852   1.64236   0.64994    0.511752    0.478862   0.00180662  0.649 ⋯
  10 │ 6.44239   1.74753   0.635937   0.615375    0.463494   0.0327256   0.635
  11 │ 8.00021   2.33972   0.689924   0.599604    1.016      0.0341892   0.689
  ⋮  │    ⋮         ⋮          ⋮          ⋮           ⋮          ⋮          ⋮  ⋱
 266 │ 4.46919   1.01337   0.16673    0.362461    0.48004    0.00414071  0.658
 267 │ 3.9543    1.00633   0.27093    0.275827    0.420701   0.0388721   0.588 ⋯
 268 │ 3.67327   0.954868  0.323962   0.282293    0.344983   0.00363011  0.633
 269 │ 2.78683   0.829576  0.394838   0.132275    0.284625   0.0178385   0.506
 270 │ 2.06009   0.643845  0.349316   0.069129    0.164816   0.0605839   0.349
 271 │ 1.87672   0.614678  0.334274   0.0691356   0.113192   0.0980763   0.334 ⋯
 272 │ 1.2605    0.437883  0.265197   0.0332794   0.0666646  0.0727419   0.265
 273 │ 1.16479   0.510662  0.236099   0.0158527   0.0708293  0.187881    0.236
 274 │ 1.31764   0.345651  0.188716   0.0119355   0.112072   0.0329278   0.188
 275 │ 1.52412   0.325618  0.155611   0.00809338  0.161913   0.0         0.155 ⋯
                                                  4 columns and 254 rows omitted, df = 275×45 DataFrame
 Row │ site     ID     IGBP     IGBPcode  IGBPname  date        year   d8     String7  Int64  String3  Int64     String3   Date        Int64  Int64 ─────┼──────────────────────────────────────────────────────────────────────────
   1 │ US-Twt     101  CRO            12  CRO       2009-01-01   2009      1   ⋯
   2 │ US-Twt     101  CRO            12  CRO       2009-01-09   2009      2
   3 │ US-Twt     101  CRO            12  CRO       2009-01-17   2009      3
   4 │ US-Twt     101  CRO            12  CRO       2009-01-25   2009      4
   5 │ US-Twt     101  CRO            12  CRO       2009-02-02   2009      5   ⋯
   6 │ US-Twt     101  CRO            12  CRO       2009-02-10   2009      6
   7 │ US-Twt     101  CRO            12  CRO       2009-02-18   2009      7
   8 │ US-Twt     101  CRO            12  CRO       2009-02-26   2009      8
   9 │ US-Twt     101  CRO            12  CRO       2009-03-06   2009      9   ⋯
  10 │ US-Twt     101  CRO            12  CRO       2009-03-14   2009     10
  11 │ US-Twt     101  CRO            12  CRO       2009-03-22   2009     11
  ⋮  │    ⋮       ⋮       ⋮        ⋮         ⋮          ⋮         ⋮      ⋮     ⋱
 266 │ US-Twt     101  CRO            12  CRO       2014-10-16   2014     37
 267 │ US-Twt     101  CRO            12  CRO       2014-10-24   2014     38   ⋯
 268 │ US-Twt     101  CRO            12  CRO       2014-11-01   2014     39
 269 │ US-Twt     101  CRO            12  CRO       2014-11-09   2014     40
 270 │ US-Twt     101  CRO            12  CRO       2014-11-17   2014     41
 271 │ US-Twt     101  CRO            12  CRO       2014-11-25   2014     42   ⋯
 272 │ US-Twt     101  CRO            12  CRO       2014-12-03   2014     43
 273 │ US-Twt     101  CRO            12  CRO       2014-12-11   2014     44
 274 │ US-Twt     101  CRO            12  CRO       2014-12-19   2014     45
 275 │ US-Twt     101  CRO            12  CRO       2014-12-27   2014     46   ⋯
                                                 37 columns and 254 rows omitted, _par = (α = 0.03265625, η = 0.069296875, g1 = 9.552734375, Am_25 = 17.671875, VPDmin = 1.21515625, VPDmax = 3.5, D0 = 0.6541015625, kQ = 0.10114375, kA = 0.89921875, S_sls = 0.01015625, fER0 = 0.152734375, hc = 0.5))
_par ==0.03265625, η=0.069296875, g1=9.552734375,
  VCmax25=17.671875, VPDmin=1.21515625, VPDmax=3.5, D0=0.6541015625,
  kQ=0.10114375, kA=0.89921875, S_sls=0.01015625, fER0=0.152734375, hc=0.5)

par = Param_PMLV2(; _par..., hc=0.5)
df.GPP_obs = df.GPPobs
df.ET_obs = df.ETobs
275-element Vector{Float64}:
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
   2.10247940297874
   3.13383220534183
   ⋮
   2.13929072217904
   1.00971020427327
   1.4844456547101
   1.60968225909663
   1.27912555152687
   1.06059787751179
   0.779535860816098
   0.578394230978662
   0.332647468807178
   0.50089104654852
   1.26067990507329
   1.35241601110533

1.1 模型参数率定

parNames = [
  :α, :η, :g1, :VCmax25, :VPDmin, :VPDmax, :D0, :kQ, :kA, :S_sls, :fER0 # :hc
]

theta, goal, flag = ModelCalib(df, par0, parNames)
df_out = PMLV2_sites(df; par=theta2par(theta, parNames))
df_out[1:10, :]
Iteration =   0, nEvals = 115, Best Cost = -0.52566
Iteration =   1, nEvals = 278, Best Cost = -0.54122
Iteration =   2, nEvals = 473, Best Cost = -0.54122
Iteration =   3, nEvals = 672, Best Cost = -0.54528
Iteration =   4, nEvals = 872, Best Cost = -0.54900
Iteration =   5, nEvals = 1073, Best Cost = -0.56046
Iteration =   6, nEvals = 1263, Best Cost = -0.56097
Iteration =   7, nEvals = 1448, Best Cost = -0.56958
Iteration =   8, nEvals = 1635, Best Cost = -0.57242
Iteration =   9, nEvals = 1832, Best Cost = -0.57700
Iteration =  10, nEvals = 2025, Best Cost = -0.57763
Iteration =  11, nEvals = 2223, Best Cost = -0.57836
Iteration =  12, nEvals = 2419, Best Cost = -0.57941
Iteration =  13, nEvals = 2616, Best Cost = -0.57968
10×14 DataFrame
Row ET GPP Ec Ecr Eca Ei Pi Es_eq Eeq ET_water Ga Gc_w fval_soil Es
Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64
1 0.414697 0.545709 0.165042 0.00883793 0.156205 0.0153008 0.313699 0.234353 0.290273 0.858553 0.014428 0.00136352 1.0 0.234353
2 0.716772 1.09093 0.541877 0.0233499 0.518527 0.0 0.0 0.265969 0.368498 1.62444 0.0162197 0.00222434 0.657576 0.174895
3 0.902411 1.11298 0.440393 0.0696786 0.370714 0.0893813 0.877869 0.372637 0.576104 1.31439 0.0106403 0.00250303 1.0 0.372637
4 1.25377 1.89176 0.80629 0.109563 0.696727 0.0168142 0.130936 0.430665 0.741149 1.8532 0.0153346 0.00400455 1.0 0.430665
5 1.51155 2.04539 0.813452 0.162363 0.651089 0.256732 1.65602 0.441369 0.84345 1.7751 0.0147179 0.00453156 1.0 0.441369
6 1.52538 1.78246 0.485462 0.103802 0.38166 0.68448 8.02602 0.355437 0.752408 1.41139 0.0268158 0.00464504 1.0 0.355437
7 1.92613 2.66826 0.889583 0.285437 0.604146 0.579148 2.99185 0.457394 1.06993 1.78514 0.0164274 0.00637722 1.0 0.457394
8 2.03736 2.9364 0.903795 0.323217 0.580577 0.712139 4.23086 0.421428 1.10398 1.75211 0.0181245 0.00727743 1.0 0.421428
9 2.39085 4.92544 1.73869 0.795952 0.942739 0.00180662 0.00669338 0.650355 1.90304 2.80294 0.0131528 0.0108613 1.0 0.650355
10 2.55767 4.9414 1.78719 0.91494 0.872252 0.134092 0.447158 0.636383 2.0679 2.85248 0.0152381 0.0113444 1.0 0.636383

1.2 拟合优度

gof = [
  (; var="ET", GOF(df.ET_obs, df_out.ET)...),
  (; var="GPP", GOF(df.GPP_obs, df_out.GPP)...)] |> DataFrame
DataFrame(gof)
2×10 DataFrame
Row var NSE R2 KGE R RMSE MAE bias bias_perc n_valid
String Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Int64
1 ET 0.615196 0.617177 0.713556 0.785606 1.34366 0.993881 -0.0778958 -2.59433 263
2 GPP 0.544157 0.620855 0.468911 0.787944 3.30045 2.62295 0.63708 16.9614 260

1.3 绘图

using Plots
gr(framestyle=:box, titlefontsize=12)
t = df.date
inds = 1:46*1

function plot_var(var; label=string(var), title=string(var),
  data=df_out, scale=1, kw...)
  plot(t[inds], data[inds, var] * scale; label, title, kw...)
end
function plot_var!(p, var; label=string(var),
  data=df_out, kw...)
  plot!(p, t[inds], data[inds, var]; label, kw...)
end

p_et = plot(title="ET components (mm/d)")
plot_var!(p_et, :Ec)
plot_var!(p_et, :Es)
plot_var!(p_et, :Ei)
plot_var!(p_et, :ETobs; data=df, label="ET_obs", color=:black)

p_gpp = plot_var(:GPP; title="GPP (gC m-2 d-1)", label="GPP")
plot_var!(p_gpp, :GPPobs; data=df, label="GPP_obs", color=:black)

plot(
  p_et, p_gpp,
  plot_var(:Eeq; title="Eeq (mm/d)", label="Eeq"),
  plot_var(:Gc_w; title="Conductance (m s-1)", label="Gc"),
  plot_var(:fval_soil; title="β_soil", label="β_soil"),
  plot_var(:VPD; data=df, scale=-1, title="-VPD (kPa)", label="-VPD"),
  plot_var(:Pi; title="P - Ei (mm/d)"),
  plot_var(:Es_eq; title="Es_eq (mm/d)"),
  size=(800, 700),
  layout=(4, 2),
)