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
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)
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 ),
)