Model background

Here is an example based on the model of Freyberg, 1988. The synthetic model is a 2-dimensional MODFLOW model with 1 layer, 40 rows, and 20 columns. The model has 2 stress periods: an initial steady-state stress period used for calibration, and a 5-year transient stress period. The calibration period uses the recharge and well flux of Freyberg, 1988; the last stress period use 25% less recharge and 25% more pumping.

The inverse problem has 761 parameters: hydraulic conductivity of each active model cell, calibration and forecast period recharge multipliers, storage and specific yield, calibration and forecast well flux for each of the six wells, and river bed conductance for each 40 cells with river-type boundary conditions. The inverse problem has 12 head obseravtions, measured at the end of the steady-state calibration period. The forecasts of interest include the sw-gw exchange flux during both stress periods (observations named sw_gw_0 and sw_gw_1), and the water level in well cell 6 located in at row 28 column 5 at the end of the stress periods (observations named or28c05_0 and or28c05_1). The forecasts are included in the Jacobian matrix as zero-weight observations. The model files, pest control file and previously-calculated jacobian matrix are in the freyberg/ folder

Freyberg, David L. "AN EXERCISE IN GROUND‐WATER MODEL CALIBRATION AND PREDICTION." Groundwater 26.3 (1988): 350-360.

The plot shows the Freyberg (1988) model domain. The colorflood is the hydraulic conductivity ($\frac{m}{d}$). Red and green cells coorespond to well-type and river-type boundary conditions. Blue dots indicate the locations of water levels used for calibration.

Using pyemu

Drawing from the prior

Now we need a prior-realized ParameterEnsemble, which stores a pandas.DataFrame under the hood.

draw

The ParameterEnsemble class has several draw type methods to generate stochastic values from (multivariate) (log) gaussian, uniform and triangular distributions. Much of what we do is predicated on the gaussian distribution, so let's use that here.

The gaussian draw accepts a cov arg which can be a pyemu.Cov instance. If this isn't passed, then the draw method constructs a diagonal covariance matrix from the parameter bounds (assuming a certain number of standard deviations represented by the distance between the bounds - the sigma_range argument)

draw also accepts a num_reals argument to specify the number of draws to make:

Note that these draw methods use initial parameter values in the control file (the Pst.parameter_data.parval1 attribute) the $\boldsymbol{\mu}$ (mean) prior parameter vector. To change that, we need to update the parameter values in the control file:

plotting

Since ParameterEnsemble stores a pandas.DataFrame, it has all the cool methods and attributes we all love. Let's compare the results of drawing from a uniform vs a gaussian distribution. The actual dataframe is stored under the private attribute ParameterEnsemble._df:

The gaussian histo go beyond the parameter bound - bad times. Luckily, ParameterEnsemble includes an enforce method to apply parameter bounds:

bayes linear monte carlo

we can use the bayes linear posterior parameter covariance matrix (aka Schur compliment) to "precondition" the realizations using linear algebra so that they hopefully yield a lower phi. The trick is we just need to pass this posterior covariance matrix to the draw method. Note this covariance matrix is the second moment of the posterior (under the FOSM assumptions) and the final parameter values is the first moment (which we parrep'ed into the Pst earlier)

Now we just need to run this preconditioned ensemble to validate the FOSM assumptions (that the realizations do yield an acceptably low phi and that the relation between parameters and forecasts is linear)