I remember looking into this problem several years ago and not finding a suitable solution. I am wondering whether FEFLOW now has a suitable method for this now?

In particular, I want to simulate nutrient mass transport in an aquifer system with the source of nutrients being diffuse, i.e. from broad application of fertilisers and leaching with irrigation and rainfall recharge. I know that in MODFLOW-MT3DMS, it is possible to set a "recharge concentration" which when combined with recharge, computes the correct mass flux to the top of the model.

With FEFLOW, I recall this had to be modelled using the mass flux (type 3) boundary. However, this could only be used in the divergent form of the transport equation and was not suitable for predicting mass discharges at type 1 constant head boundaries. Using the convergent form restricts the user to define type 1 mass flux boundaries, which are not able to simulate the desired mass flux.

Is there a practical solution to this problem using FEFLOW?

Cheers,
Blair
This is a common question in mass transport problem in FEFLOW. Unfortunately, there is no a single approach to answer, instead the conceptualization of the mass-transport boundaries and mass-transport formulation (convective/divergence) in FEFLOW should be addressed at the project scale.
Below I described a non-common boundary configuration, which could help you. If the task is just to apply a “recharge concentration”, this means a combined specification of water inflow and mass load (dispersive and advective) at the same location. The boundary conditions must ensure that the total of dispersive and advective flux equals the prescribed mass load at all times (i.e. your conceptual amount). The best combination is a Mass-Transfer BC (i.e. Cauchy-type and equal to the desired concentration) and a Fluid-Flux BC (Neumann-type) with the default mass transport formulation. The transfer-rate for the Cauchy-BC should be equal to the incoming fluid flux. The fact you are using the convective form of transport allows you a free outflow from the model at outflowing borders. This boundary configuration is valid for weak/strong lateral hydraulic gradients.
Since you mentioned fertilization problems, I would like to mention we have recently extended the FEFLOW capabilities for reactive transport with our FEFLOW plug-in piChem, which is able to couple FEFLOW models (2D/3D) with the reactive engine PhreeqC.
Hi, I am recently using piChem to simulate rock water reactiions, I have a very complicated structure with 5 aquifers and 4 aquitards.  I tried to define each zone as single reactive element, but after I define all these zones attached to phreeqc files respectively, feflow seems to have crash problem, the time run panel buttons are unacitve all the time, and the model does not run at all, always saying 'computing' but the time step stays at 0d without any other errors. The worst thing is, I can not even close feflow ! But when I tried the piChem examples, or establish a simple model with same definitions using phreeqc, for example, one layer model, it could run. Therefore, does anybody know how I should solve it? Is it because the model layers are too many for the computer to calculate or? any suggestions to improve?