Hi Leelja,
I would suggest to put the 35000 mg/L as initial only at the seaside part of the domain and have 0 mg/L in the rest of the domain, likely where you know there's only freshwater.
Also, all inflowing BC's on the freshwater side should have a 0 mg/L mass BC.
Apart from the mass transport, I would also suggest to run the model as flow only using the van Genuchten, just to make sure this is stable. If the head solution is already unstable, you have to fix that first before working on the transport.
When the model is being slow with constraints, I would assume they are frequently being triggered during the run. You can easily check that by saving a .dac file and plot all mass BC's. Then you can cycle through each time-step and check. A frequent - and maybe not desired - triggering of constraints can often be due to instability of the model, i.e. over- or undershooting of the primary solution (i.e. mass concentration). If this is the case, you should try stabilizing the model by adjusting the mesh first, then work on other parameters like dispersivities or solver settings. A general good idea with such a combinded unsaturated and mass transport approach is working with the FE/BE time-stepping scheme. Further, try putting a growth factor (starting at small values 1.05 to 2, you have to do test runs for optimal behaviour).
Since it is being unstable when using van Genuchten, I would also suggest to have a look at the saturation and its distribution in each time-step. You might also have to work on the mesh in order to better represent the saturation, or choose different van Genuchten parameters for adjusting the saturation solution.
