def cf_model(model, k, i, j, base, Q=-100):
    wd = {1: [[k, i, j, Q]]}
    model.remove_package('WEL')
    wel = flopy.modflow.ModflowWel(model=model, stress_period_data=wd)
    wel.write_file()
    model.run_model(silent=True)
    # get the results
    hedObj = flopy.utils.HeadFile(os.path.join(cf_pth, 'DG.hds'), precision='double')
    cbcObj = flopy.utils.CellBudgetFile(os.path.join(cf_pth, 'DG.cbc'), precision='double')
    kk = hedObj.get_kstpkper()
    v = np.zeros((len(kk)), dtype=np.float)
    h = hedObj.get_ts((k, i, j))
    for idx, kon in enumerate(kk):
        if h[idx, 1] == model.lpf.hdry:
            v[idx] = np.nan
        else:
            v1 = cbcObj.get_data(kstpkper=kon, text='DRAINS', full3D=True)[0]
            v2 = cbcObj.get_data(kstpkper=kon, text='STREAM LEAKAGE', full3D=True)[0]
            v3 = cbcObj.get_data(kstpkper=kon, text='ET', full3D=True)[0]
            v[idx] = ((v1.sum() + v2.sum() + v3.sum()) - base) / (-Q)
    return v
dataset/ETHPy150Open modflowpy/flopy/examples/groundwater_paper/scripts/uspb_capture.py/cf_model/original