Possible bug in RCT

[emorway-usgs]

An email from Paul Hsieh:
"The problem seems to occur when sorption is involved in a multi-layer model. I created some runs to illustrate this. The set up consists of 4 layers, 1 row, 10 columns, with flow occurring along the row from left to right. The set ups were created with ModelMuse.

When there is no sorption, both mt3d-usgs and mt3dms give the same result.

However, when sorption (linear isotherm) is turned on, mt3d-usgs does not appear to correctly read the sorption parameters (bulk density and Kd), so the retardation factor is not correctly computed. (See the listing file example.mls in the section "SORPTION AND 1ST/0TH ORDER REACTION PARAMETERS". The retardation factor should be 1.8 in all layers.)

Running the same sorption problem using mt3dms appears to give the correct result."

[emorway-usgs]

Richard Winston had a look in ModelMuse and the source code, and responded: "I think I found the issue. Data set 1 of the reaction package in MT3D-USGS includes an extra variable, IREACTION, that is not present in in MT3DMS. ModelMuse does not write a value for that variable. Because that variable is not present, MT3D-USGS activates the following code starting on line 51 of rct1.f.
100 IF(IERR.NE.0) THEN
IRCTOP=1
IGETSC=0
BACKSPACE (INRCT)
READ(INRCT,'(2I10)') ISOTHM,IREACT
ENDIF
This code sets IRCTOP to 1 which is different from the value specified in the input file.

This same code is present in MT3DMS but doesn't get activated there because MT3DMS does not attempt to read IREACTION. I suspect this code is present for backwards compatibility with MT3D.

I think the best way to handle this would be to insert a similar if statement to the one above that would read only 4 variables instead of 5. It should supply a default value for IREACTION. If an error occurs in reading the 4 variables, the above if statement can be used for backwards compatibility with MT3D.

I will modify ModelMuse to write a default value for IREACTION until I can get to more fully supporting MT3D-USGS."

[emorway-usgs]

A solution was offered by Richard Winston and adopted with commit b15eaa4

Here is Richard's suggestion:

I have updated ModelMuse to write a default value of 0 for IREACTION. you can get the updated version from the following URLs.
32-bit version: ftp://ftpext.usgs.gov/pub/er/va/reston/rbwinst/ModelMuse32_3_9_0_24.zip
64-bit version: ftp://ftpext.usgs.gov/pub/er/va/reston/rbwinst/ModelMuse64_3_9_0_24.zip
To install, replace your existing version of ModelMuse.exe with the one from the zip file.

The following is my suggested change to rct1.f starting at line 48.

C--READ AND ECHO SORPTION ISOTHERM TYPE AND FLAG IREACT
READ(INRCT,'(6I10)',ERR=100,IOSTAT=IERR)
& ISOTHM,IREACT,IRCTOP,IGETSC,IREACTION
100 IF(IERR.NE.0) THEN
! Backwards compatibility with MT3DMS.
IREACTION=0
IERR = 0
BACKSPACE (INRCT)
READ(INRCT,'(4I10)',ERR=101,IOSTAT=IERR)
& ISOTHM,IREACT,IRCTOP,IGETSC
ENDIF
101 IF(IERR.NE.0) THEN
IRCTOP=1
IGETSC=0
IREACTION=0
BACKSPACE (INRCT)
READ(INRCT,'(2I10)') ISOTHM,IREACT
ENDIF