Galilee Basin hydrogeological model review

A regional numerical groundwater model of the Galilee Basin – here referred to as the Galilee Basin hydrogeological (GBH) model – has been developed utilising data and interpretation undertaken for the Galilee subregion by the Assessment team. The model was jointly funded by the Commonwealth Department of the Environment’s Office of Water Science and the Queensland Government. A consultant – HydroSimulations – was engaged to produce a calibrated and stress-tested transient numerical groundwater flow model capable of simulating the cumulative impacts of proposed coal mining developments. The GBH model was built over a 2.5 month period. The following summaries are based on a detailed description of the GBH model, provided in Turvey et al. (2015).

Input data for the model was compiled as part of this assessment for the Galilee subregion. Companion product 2.1-2.2 for the Galilee subregion (Evans et al., 2018a) provides detail on data used as input to the GBH. These data include: the geological model of the Galilee subregion, potentiometric surfaces for aquifers in the subregion, and recharge volumes, among other datasets. The GBH model was built utilising the new MODFLOW-USG platform and AlgoMesh, which uses flexible meshes to allow refinement around key features while maintaining a reasonable cell count for computational efficiency. The GBH model consists of an unstructured grid (except for mines) consisting of 1.93 million cells ranging in size from 150 m to 10 km and covers an area of 299,400 km2. A cell size of 165 m to 250 m is used around modelled mines. Modelling was carried out in three stages, resulting in steady-state, transient and predictive simulations. Mines simulated in the transient model are: China First, Kevin’s Corner, Carmichael, South Galilee, China Stone, Hyde Park and Alpha. These are mines that comprise the modelled CRDP, as outlined in companion product 2.3 for the Galilee subregion (Evans et al., 2018b).

The GBH model consists of 13 layers representing the major hydrostratigraphic units present within the subregion (Table 6). In many cases, individual model layers represent multiple stratigraphic units that are grouped based on similar hydraulic properties.

Table 6 Model layer thickness and corresponding hydrostratigraphic groupings (after Turvey et al., 2015)

Model layer

Hydrostratigraphic unit


Mean thickness


Maximum thickness


Calibrated horizontal hydraulic conductivity (Kh) (m/day) (same for steady-state and transient calibration)

Specific storage (Ss)

(transient calibration)


Specific yield (Sy)

(transient calibration)


Layer 1


Unconfined aquifer



1.00x100 (regolith) /3.00x101 (Cenozoic)


1.00x10–1 (Cenozoic) /1.00x10–2 (regolith)

Layer 2

Winton-Mackunda formations

Unconfined to semi-confined aquifer






Layer 3

Allaru/Toolebuc/ Wallumbilla formations (Rolling Downs Group aquitard)







Layer 4

Cadna-owie Formation/Hooray Sandstone/Gilbert River Formation/ Ronlow beds

Confined aquifer





1.00x10–2 (east) /8.00x10-3 (west)

Layer 5

Westbourne/Adori/ Birkhead formations (Injune Creek Group)







Layer 6

Hutton Sandstone/Boxvale Sandstone /Precipice Sandstone

Confined aquifer



1.27x10–1 (east) /7.59x10–1 (west)


1.50x10–2 (east) /7.00x10–3 (west)

Layer 7

Moolayember Formation







Layer 8

Warang-Clematis Group

Confined aquifer



1.96x100 (east) /4.02x100 (west)


1.50x10–2 (east) /6.00x10–3 (west)

Layer 9

Rewan Group (includes Dunda beds)







Layer 10

Top Bandanna Formation to base of B Seam (BC1)

Partial aquifer






Layer 11

Base B Seam to Top of E seam (BC2)

Partial aquifer






Layer 12

Top of E seam to base of Colinlea Sandstone (BC3)

Partial aquifer






Layer 13

Joe Joe Group/Drummond Basin sediments/crystalline basement




8.77x10–5 (Joe Joe Group) /2.00x10–3 (Drummond Basin) /2.00x10–5 (crystalline basement)



Mean thickness and maximum thickness figures are from Table 3-1 in Turvey et al. (2015).

Significant vertical hydraulic gradients exist through the upper Permian coal measures. For this reason the upper Permian coal measures were modelled as three layers (layers 10, 11 and 12). From top to bottom these sub-units are informally designated BC1, BC2 and BC3. The partitioning was done on the basis of approximately similar groundwater pressures existing in each sub-unit (Table 7). The upper Permian coal measures are designated as a partial aquifer because of low hydraulic conductivity (variable Kh values, but less than 0.2 m/day).

Table 7 Stratigraphic units and marker beds used to define subdivision of the upper Permian coal measures



Marker bed


Bandanna Formation / Blackwater Group

Top Bandanna Formation to base B seam


Peawaddy Formation / Black Alley Shale / Colinlea Sandstone (upper)

Base B seam to base DE Sandstone


Colinlea Sandstone (lower)

Top E seam to base Colinlea

Major rivers within the subregion are represented in the model using the MODFLOW River package and are assumed to be connected to the groundwater system. Rivers modelled include: Carmichael River, Belyando River, Native Companion Creek, Alpha Creek, Aramac Creek, Thomson River, Flinders River, Diamantina River, Barcoo River, Warrego River, and Blackwater Creek.

Groundwater discharges from 121 springs are simulated in the model using the MODFLOW River package (head-dependent boundary conditions) by setting the stage and bed elevation of the river cells to be equal, allowing one-way flow (discharge from groundwater) only. The River package allows for the specification of spring – discharge observation points where the fluxes from these river cells are functions of the calculated groundwater level and assigned conductance. Major springs represented include Doongmabulla, Mellaluka, Edgbaston, Coreena and Corinda springs. Based on limited data available from the Queensland springs database, the major springs were assigned flows greater than 100 m3/day; all other springs modelled were limited to less than 100 m3/day.

Long-term mean recharge rates used in the model are based on chloride mass balance estimates reported in companion product 2.1-2.2 for the Galilee subregion (Evans et al., 2018a). Estimated recharge rates vary from a low of 0.1 mm/year, for the Joe Joe Group, to a high of 5 mm/year for Cenozoic alluvial deposits.

It is considered that no lateral groundwater inflow or outflow occurs along the eastern margin of the Galilee subregion. Consequently, no boundary condition is assigned to the eastern edge of all model layers. In contrast, general-head boundaries are applied along the western boundary of the model to the aquifer model layers 2, 4 and 6, representing the regional hydraulic gradient and communication with aquifers outside of the model domain.

Existing groundwater extraction is represented in the model using the MODFLOW Well package. For steady-state simulations the mean recorded artesian well discharge between 1983 and 2012 was applied to be consistent with the data used for calibration. For the historical transient model the actual annual discharge rate values were used, while predictive model runs used a constant flow rate of 895 ML/day, which represents estimated total daily extraction for the year 2012.

Turvey et al. (2015) noted that bore discharge estimates provided may be inaccurate and that for some bores there is uncertainty about the source aquifer. As a result it is possible that well stresses are misapplied in the model. To compensate, variable bore discharge rates have been applied during model simulation to prevent incorrectly assigned discharge rates from drying out cells in low conductivity layers, causing model instability.

Mine dewatering for the seven mines included in the predictive model is simulated using the MODFLOW Drain package. The package integrates the proposed mine plan location, timing, depth and method of extraction. Open-cut mine pits are assumed to be the full thickness of the relevant target coal seams. Drain cells are progressively added to the transient model to simulate extraction of mine panels over time.

Longwall mine induced fracturing is simulated using the MODFLOW-USG Time-Variant Materials package. Hydraulic parameters are changed with time in the goaf and overlying fracture zones during extraction of each longwall panel. Steady-state calibration

Steady-state calibration was assessed against water levels corrected for temperature and salinity to ensure comparability with MODFLOW predictions which do not take account of density variation.

Calibration statistics for the steady-state simulation are shown in Table 8. Predictions within ±10 m of target levels were distributed evenly across the model domain, however a tendency remains to under predict water levels in the west and over predict water levels in the east. An area of significant over prediction (>50 m) occurs within layers 11 and 12 at the south-east of the model. Target residuals at the southern group of mines (Kevin’s Corner, Alpha, China First and South Galilee) are relatively good, with an overprediction of approximately 13 m being the most significant discrepancy. However, at Carmichael, particularly north of the Carmichael River, significant model overpredictions of up to 38 m occur within the coal seam layers. No data were available for calibration at either China Stone or Hyde Park mines at the time of model development.

Table 8 Steady-state calibration statistics (after Turvey et al., 2015)



Residual mean (m)


Absolute residual mean (m)


Residual standard deviation (m)


Sum of squares (m2)


Root mean square (RMS) error (m)


Minimum residual (m)


Maximum residual (m)


Number of observations


Range in observations (m)


Scaled RMS error


% Targets within ±10 m


% Targets within ±25 m


Turvey et al. (2015) includes a comparison of simulated potentiometric surface contours derived from the steady-state calibration and contours interpreted from observed water levels for companion product 2.1-2.2 for the Galilee subregion (Evans et al., 2018a). Turvey et al. (2015) found that there is general agreement between the pattern of the Eromanga Basin modelled surfaces and interpreted surfaces. A comparison of model layers 1 and 2 and interpreted surfaces for the Cenozoic and Winton and Mackunda formations shows differences in elevation generally less than 20 m and similar inferred flow directions. A comparison between the model layer 4 surface and the interpreted surface for the Hooray Sandstone showed modelled elevations 10 to 30 m below the interpreted surface and similar inferred flow directions. These properties are the same for the comparison between the model layer 6 surface and the interpreted surface for the Hutton Sandstone.

The discrepancies between the modelled Galilee Basin units (model layers 8, 10, 11 and 12) and interpreted surfaces are more variable with a general pattern of under prediction of modelled water levels over most of the model domain and over prediction of modelled water levels toward the eastern boundary of the subregion where aquifers outcrop. However, inferred flow directions are generally similar for both modelled and interpreted surfaces.

The general underestimation of modelled water levels in the location of some springs is likely to account for predicted spring flows that are less than observed flows. The total modelled spring flow under estimates reported flow by around one order of magnitude.

Calibration to baseflow was not undertaken due to uncertainty in available baseflow records.

The water balance for the steady-state simulation indicates that just under half of the recharge to groundwater is from direct rainfall recharge, with an approximately equal contribution from river leakage and inflow along model boundaries accounting for the remaining portion of inflow. Losses from the system in order of magnitude are: evapotranspiration; groundwater extraction; discharge to rivers and springs; followed by groundwater outflow along model boundaries. All but two rivers in the model (Thomson and Diamantina rivers) act as net losing systems overall. Transient calibration

Manual transient calibration was carried out for the period January 1998 to December 2012. Due to the limited data available for transient calibration, steady-state Kh and Kv values were not varied and only variation in storage parameters (Ss and Sy) was allowed. High specific storage (Ss) values were required to prevent excessive drawdown due to well extraction over the calibration period. For the transient simulation, groundwater recharge rates were varied using a series of recharge factors for the period 1983 to 2011 derived from the Australian Water Resources Assessment landscape model (AWRA-L) in order to scale steady-state recharge values to derive a transient recharge dataset. Recharge estimates for the period 2012 to 2014 were derived using the Penman-Grindley recharge model (Finch, 1994).

Calibration statistics for the transient simulation are shown in Table 9.

Table 9 Transient model calibration statistics (after Turvey et al., 2015)



Residual mean (m)


Absolute residual mean (m)


Residual standard deviation (m)


Sum of squares (m2)


Root mean square (RMS) error (m)


Minimum residual (m)


Maximum residual (m)


Number of observations


Range in observations (m)


Scaled RMS error


% Targets within ±10 m


% Targets within ±25 m


The transient model mass balance indicates that the majority of water that enters the model does so via a combination of river leakage and rainfall recharge. Regional groundwater flow into the model, via general-head boundaries, is approximately double that of outflow. Over the calibration period there was a net loss in storage suggesting a depleting resource. Predictive runs

Two predictive model runs were simulated and run out to the year 2220 in order to allow for sufficient simulation of the water level recovery phase.

The first predictive model run for the GBH model represents a baseline which consists of existing groundwater use without groundwater extraction due to coal mining (Turvey et al., 2015). For BA, the baseline is defined as a future that includes all coal mines and coal seam gas (CSG) fields that are commercially producing as of December 2012 (see Section 2.3.4 in companion product 2.3 (Evans et al., 2018b) for further detail). As of December 2012, there were no operational coal mines or CSG fields in the Galilee subregion (Evans et al., 2018b). Hence, the first predictive run of the GBH model represents the coal resource development baseline defined for this assessment.

The second predictive model run includes the seven coal mines that comprise the modelled CRDP (see Section in companion product 2.3 (Evans et al., 2018b) for detail) and represents the cumulative drawdown of these potential coal mining operations. Additional stresses to those imposed on the baseline predictive model run include: a progressive excavation and associated dewatering of the seven open-cut and longwall mines; progressive development of fracture zones above underground mines and associated changes in permeability; and progressive placement of fill in open-cut mines to end of mine life using modified hydraulic parameters and recharge. Limitations – Galilee Basin hydrogeological model

The GBH model is a complex regional groundwater model produced in a very short time period and within a limited budget. Upon completion, a number of limitations and issues have been identified by the contractors – HydroSimulations – and reviewers from Queensland Government, Office of Water Science and relevant project and discipline leads from the Bioregional Assessment Technical Programme. The limitations noted by various reviewers include:

  • The ability of the model to simulate impacts to the unconfined Cenozoic aquifer – model layer 1 – is limited. Predictive model runs did not show any impact to model layer 1, particularly in the vicinity of open-cut mine pits. It would be expected that drawdown and dewatering of Cenozoic units would occur adjacent to mine pits which penetrate and expose the Cenozoic units. Turvey et al. (2015) suggested that certain options within the MODFLOW-USG code do not allow the model to represent negative pressures (i.e. dry cells representing unsaturated conditions). This means that an artificial head is maintained at the base of each layer, which may result in under prediction of the magnitude of drawdown within each cell. This will have had particular implication for impacts predicted in model layer 1, essentially resulting in water being held at the base of model layer 1 for the entire model simulation and resulting in negligible reported drawdown within this layer.
  • Combining several hydrogeological units into a single model layer means that accurate representation of real-world groundwater levels for all members of the GAB and Galilee Basin units cannot be produced. This is a particular issue for the upper Permian coal measure layers (model layers 10, 11 and 12), where bulk hydraulic conductivity parameters applied to the model representing coal seam layers are higher than those of the interburden layers (that are not included). Additionally, no vertical hydraulic conductivity variation exists within the bulked layers. As discussed in Section, combining hydrogeological layers together can be justified when conservatively estimating drawdowns. The lumping of hydrogeological layers, however, can result in biased parameter estimates if local groundwater level observations are used to constrain the aggregated hydraulic parameters. Another potential issue is that mine dewatering rates will be over-estimated as the aggregated unit will be dewatered in the model, rather than just the mined interval. Such an overestimate can be considered conservative and therefore justified. It will contribute to discrepancies between estimated pumping rates between models, further complicating comparisons between model results.
  • The representation of losing and gaining streams in the GBH model is suited to predict the hydrological change due to coal resource development. However, in a deterministic aquifer simulator such as the GBH model, incorrectly specified surface water – groundwater boundary conditions is a large source of conceptual model uncertainty. This may bias parameter estimates and therefore bias predictions. Debate exists about the representation of the Thomson, Alice, Carmichael and Belyando rivers as losing streams. Initial assessments done for companion product 2.1-2.2 for the Galilee subregion (Evans et al., 2018a) suggests these streams could be predominantly gaining. Section identifies the large potential of improving the characterisation of the connection status of the river network to reduce predictive uncertainty.
  • Queensland Government reviewers noted issues with some of the input data, in particular the boundary definitions for some model layers (e.g. the Rewan Group) and the accuracy of some well flow rates. There is a possibility that over extraction of groundwater from the model due to erroneous input well flow rate data has resulted in calibrated hydraulic conductivities and storage properties being too high. Anomalously high hydraulic conductivities and storage properties are likely to result in increased mine inflows and greater resultant drawdowns.
  • Limited data were available for calibration; most data collected were of single-point observations over a period of 30 years. Therefore, reliability of the data is unknown as long-term trends were not able to be analysed. Section highlights the need for additional field data collection.
  • Limited formal sensitivity analysis has been carried out due to the short model development time frame and agreed scope for the project. It is therefore difficult to clearly indicate which parameters are controlling calibration, inflows and predicted impacts.
  • The GBH model input data does not include more recent groundwater data that became available after the release of the groundwater modelling report for China Stone Coal Project. Extra data in AGE Pty Ltd (2015) would refine the results of the steady-state and transient models in the vicinity of China Stone and Carmichael coal projects.

Last updated:
17 December 2018
Thumbnail of the Galilee subregion

Product Finalisation date