Simulation Log

1) First simulation(s) are to establish working model and get some idea of scale of values needed for the recharge, hydraulic conductivity, and saturated thickness.  

The values that resulted in stable (no dry cells, convergence to a solution in all time steps) are:

Ksat = 3.0 meters/year (seems quite low)
Sy   = 0.5  (a bit large for geologic materials; probably 0.35 would be more defendable)
Start Head = 1822.00 meters
Bottom Elevation = 1500.00 meters
Recharge = 0.02 meters/year (20 mm/year)
Constant Head Boundary Value = 1822.00 meters
All units are meters and years.  All wells are inactive.

2) Second simulation, replace the constant recharge with a distributed value that is computed as:  

0.00004 * (Annual Rainfall - Annual Evapotranspiration) on a cell-by-cell basis.  The constant scales the values from millimeters per year into meters per year, as well as adjusts the values so the maximum values do not exceede 20mm/year and clobber the model.

3) Third simulation, activate the 9 wells.

Simulation is stable.  Ready to adjust the output control to get the equilibrium heads before wells and enabled.
Construced R script to kind of read the MODFLOW output file and make some graphics.  

Next step is to make a script to convert well locations to grid coordinates (COL,ROW)

Reading the output files is annoying, may need to modify the MODFLOW code to write a formatted structure or
write a Python script to read the output and structure for plotting.  

The Python solution is probably more maintainable in the long run.