SetConcentrations correctly in 2D model

[Ellis]

Hello David,

I am trying to model reactions in a 2D domain using PhreeqcRM. My 1D transport model works fine with PhreeqcRM, while my 2D model works without PhreeqcRM. However, I get an “unknown exception” error (I am using PhreeqcMatlab – with the advection_cpp file as my guide) when I couple the 2D model with PhreeqcRM. I know there may likely be a syntax problem going on, but I can’t figure out how to fix it. That is, I am not sure in what form I need to set the concentration before Runcells().
•   For instance, in 1D, the initial concentration at time_step 0 is an n x ncomps size array, where n = number of cells, and ncomps = number of components. Should this remain the same in 2D, with n = total number of 2D cells? Should this be the same size for the rest of the time_steps?
•   What about the mapping? I know grid2chem is 0:nxyz-1 in 1D. Not sure how it’ll look like in 2D - or if it really does matter?
•   What about ic1, ic2, f1, bc1, bc2, bc_f1, bc_conc? Are their array sizes expected to change from what they are in 1D?

In a nut shell, what is the key to setting the 2D concentration (or how should the set concentration look like) before employing Runcells()? This is where my code fails.

Thanks for any leads.

[dlparkhurst]

For starters, when you create an instance of PhreeqcRM, you must define nxyz, which is the total number of cells in your transport model. If it is a rectangular grid, nxyz = nrows * ncolumns.

You then must decide on a numbering system, 1 to nxyz, for the transport grid. For Fortran, the numbering is usually 1,1 (1); 2,1 (2); 3,1 (3), ... . For C/C++ it is usually 1,1 (1); 1,2 (2); 1,3 (3) ... . I don't know the order of 2D arrays in Matlab. You will have to make sure the component concentrations  end up in the right order to be able to transport each component. Whatever you pick, by default, there are nxyz chemistry cells with the same numbering as the transport cells. If there are inactive cells, or symmetry, you can use a mapping to limit the number of cells for which chemistry is run (nchem <= nxyz).

The arguments to SetConcentrations and GetConcentrations are always arrays that are nxyz x ncomps, regardless of whether or not there are inactive cells. The concentration data is translated from the nxyz transport cells to the nchem chemistry cells (SetConcentrations) and from the chemistry cells back to the transport cells (GetConcetrations).

Similarly, ic1, ic2, ... arrays are nxyz x 7. 

[Zhaoyang]

Hi David,

I have a followed-up question about defining the boundary cell in the PhreeqcRM input file. For 1D simulations, we usually define a solution 0 as the boundary cell. How should we do for 2D cases with more boundary cells? Thanks for your help in advance.

Best regards,
Zhaoyang

[dlparkhurst]

It just depends on the physical situation that you are modeling. If you assume uniform flow and chemical conditions in a 2-D system, you would apply solution 0 all along the boundary. If you are modeling a reactive wall, where flow is restricted to a central section of the boundary, then you would define a boundary condition only on the central section.

You could also conceive of more complicated boundaries with flow or chemistry varying across the length of the boundary. The boundary conditions may be constant head, specified flux, or head-dependent leaky. In most cases, the heads or flux associated with the boundary are time-wise constant or linearly varying in a specified direction. Similarly, the chemical condition is usually  defined as time-wise constant or linearly varying between two compositions in a specified coordinate direction.