Macropore flow solution method

The numerical formulation has to take into account the fact that flows in the macropore and in the matrix domains occur with significantly different veloci­ties. Priority in development of the numerical scheme has been put on pre­serving the water balance and ensuring numerical stability at time steps that are not much lower than the time steps used in solving the matrix flow equa­tion (Richards) alone. The solution method is mass conserving. The time step length is controlled by specifying certain limits for flow and exchange flow (depths) per time step, partly for ensuring correct dynamics of the macropore flow description and partly to avoid instability of the Richards solution due to high source/sink terms. The time step is controlled by performing an extra (a priori) macropore computation at the start of each matrix flow time step—with the matrix conditions from the previous time step. In case the resulting maxi­mum (a priori) macropore flows, infiltrations, and exchange flows exceed the specified limits, a reduced time step is estimated, assuming linear relation­ship between time step length and flow volumes (unchanged flow rates). The procedure is repeated until the estimated maximum flows are within limits.

After the time step check, a normal time step simulation is performed, solving the Richards equation for the matrix flow—with reduced source/ sink terms from macropore-matrix exchange flows of the previous time step. After this the corresponding macropore time step is performed.

The calculation procedure consists of a double sweep algorithm with the fol­lowing characteristics:

Downwards sweep (sweep 1)

First, a downwards sweep is performed, i.e. downwards flow from each cell to the cell below (or ground water table) in the macropore and in the matrix domains, and exchange flows, are calculated. Mass conservation is ensured by reducing the outflow from a cell if it exceeds the storage volume of that particular cell. The flow is not influenced by the water content of the receiving cell below, but the flow is set to 0 if the cell below has no macropores.

In the downwards sweep the downwards flow from a cell is calculated as the average of two estimates, and exchange flow as an average of four esti­mates:

·         First, flow and exchange flow estimates are calculated as function of the macropore water content at the start of the time step. The estimates are reduced, if they exceed the start volume plus the inflow from the cell above.

·         The second estimate is the flow and exchange flows as a function of the resulting macropore water content from the first estimate, including the inflow from the cell above (or the macropore infiltration if it is the upper cell). Again a reduction is made if the resulting volume would become negative.

·         The final flow and exchange flow estimates are calculated as the aver­age of the first two estimates. The flow is used as inflow to the cell below (or macropore recharge to groundwater if it is the lowest cell).

As mentioned above, each of the two macropore-matrix exchange flow esti­mates are calculated as two sub-estimates: The first estimate is calculated as function of the matrix water content at the start of the time step, and the sec­ond as function of the matrix water content as result of the first estimate. The same macropore water content (pressure) is used for both exchange esti­mates. Each sub-estimate is limited by two conditions:

·         Flow from macropore to matrix is reduced if the resulting matrix water content would exceed saturation; and

·         Flow from matrix to macropore is reduced if the resulting matrix pressure would be below the macropore pressure.

The macropore pressure used for exchange calculation is calculated as the macropore saturation (i.e. actual water content divided by porosity) multiplied by the cell height, and then reduced by the entry pressure. If the macropores of the cell are fully saturated, the macropore pressure of the cell above is added (hydrostatic conditions).

Upwards sweep (sweep 2)

In the second, upwards sweep, the macropore flows and the matrix-macrop­ore exchange flows are reduced for situations where the macropores would otherwise become over-saturated. The resulting macropore water contents from downward sweep are checked, and when they exceed the macropore porosity the flow from the cell above and the exchange inflow from the matrix are reduced accordingly. If the cell receives exchange inflow from the matrix, the exchange flow and the flow from the cell above are reduced by the same proportion, otherwise only the flow from the cell above is reduced. Mass con­servation is ensured by adding the flow reduction volume to the volume of the cell above (converted to water content).

The calculational time step is automatically adjusted in situations with high macropore and exchange flows so that these flows do not exceed certain lim­its, which from experience is known to generate numerical instability. In case of an increase of the ground water table, the water content in macropores now located below the groundwater table will be added to the groundwater recharge.

Notes

·         Macropore porosity is added to the saturated water content to get the total porosity of the cell. This is important when adjusting the SZ-UZ exchange and the calculation of the relevant specific yield correction.

·         The retention curve is only use for the matrix flow. It is not used for the macropore flow.