Matrix permeability is a key factor in determining long term gas production from shale reservoirs - requiring that it is determined under true reservoir conditions. We suggest a variable pressure gradient (VPG) protocol to measure shale matrix permeability using real reservoir fluids in powdered samples. The VPG method is described and a mathematical protocol for its analysis is developed. The first measures gas fractional production rate history under constant external pressure for each production stage and with a designated pressure gradient. The second establishes the mathematical protocol for analysis using pseudo-pressure to accommodate both the effect of gas pressure-dependent PVT parameters and desorption rate coefficient. The matrix permeability is determined by matching the solution of the model with the experimental data. The model fits the experimental data well when the fractional production is < 0.75. Shale matrix permeability is calculated in the order of magnitude of 10(-7)-10(-6) md. Methane permeability decreases with a decrease in both average pore pressure and particle size of the individual component grains. Permeability considerably more sensitive to changes in desorption rate coefficient than flow regimes. Compared with current small pressure gradient (SPG) methods, the VPG method is considerably more applicable to actual gas production and reduces to the SPG method under simplified boundary conditions. Although some approximate treatments are used for establishing the VPG method and some flow mechanisms are not considered, this study still provides an information-rich technique to determine shale matrix permeability at conditions close to reality.