I hope it helps, I am new to using finite element mesh. Is there a source I could look at for the diffusion equation you are modeling? Probably there is a better way to model a single source in Paris than the Gaussian used in the second plot.

I don't know if this solution is what you are going for - it is more of a kludge, and I can't say if it correctly models diffusion from a center to a regions boundaries.

Method 1 - add an interior boundary and second DirichletCondition

But this is a quick method. First you use the "IncludePoints" option of ToBoundaryMesh to include a set of points describing a circular region around Paris. These points are then considered part of the boundary for NDSolve. Now you add a second DirichletCondition setting the value of u to some high constant at that interior boundary.

The two adjustable parameters here are the radius around Paris that you consider, and the value you set the function to at that boundary.

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
parisradius = 0.3;
parisdiffusionvalue = 150;
carto = DiscretizeGraphics[
   CountryData["France", {"Polygon", "Mercator"}]];
paris = First@
   GeoGridPosition[
    GeoPosition[
     CityData[Entity["City", {"Paris", "IleDeFrance", "France"}]][
      "Coordinates"]], "Mercator"];
bmesh = ToBoundaryMesh[carto, AccuracyGoal -> 1, 
   "IncludePoints" -> CirclePoints[paris, parisradius, 50]];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 5];
mesh["Wireframe"];
op = -Laplacian[u[x, y], {x, y}] - 20;
usol = NDSolveValue[{op == 1,
    DirichletCondition[u[x, y] == 0, Norm[{x, y} - paris] > .6], 
    DirichletCondition[u[x, y] == parisdiffusionvalue, 
     Norm[{x, y} - paris] < .6]}, u, {x, y} \[Element] mesh];
Show[ContourPlot[usol[x, y], {x, y} \[Element] mesh, 
   PlotPoints -> 100, ColorFunction -> "Temperature"], 
  bmesh["Wireframe"]] // Quiet
Plot3D[usol[x, y], {x, y} \[Element] mesh, PlotTheme -> "Detailed", 
 ColorFunction -> "Rainbow", PlotPoints -> 50]

Thanks to kglr for showing how to convert city coordinates to the Mercator projection.

Method 2 - Modify the differential equations

This seems to give a better result, but again it isn't based on any transport theory. All I'm doing is setting the Laplacian equal to a Gaussian function centered around Paris. The only important adjustable parameter is the width of the Gaussian,

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
parisradius = 2;
parisdiffusionvalue = 150;
carto = DiscretizeGraphics[
   CountryData["France", {"Polygon", "Mercator"}]];
paris = First@
   GeoGridPosition[
    GeoPosition[
     CityData[Entity["City", {"Paris", "IleDeFrance", "France"}]][
      "Coordinates"]], "Mercator"];
bmesh = ToBoundaryMesh[carto, AccuracyGoal -> 1];
mesh = ToElementMesh[bmesh, MaxCellMeasure -> 5];
mesh["Wireframe"];
op = -Laplacian[u[x, y], {x, y}];
usol = NDSolveValue[{op == 
     150 Exp[-Norm[{x, y} - paris]^2/parisradius],
    DirichletCondition[u[x, y] == 0, True]}, 
   u, {x, y} \[Element] mesh];
Show[ContourPlot[usol[x, y], {x, y} \[Element] mesh, 
   PlotPoints -> 100, ColorFunction -> "Temperature"], 
  bmesh["Wireframe"]] // Quiet
Plot3D[usol[x, y], {x, y} \[Element] mesh, PlotTheme -> "Detailed", 
 ColorFunction -> "Rainbow", PlotPoints -> 50]

Ali - I don't know what happened here, I was most certainly using version 10.xx - not sure which, but when I evaluate this now it doesn't work. The problem is the one mentioned by the OP here: http://community.wolfram.com/groups/-/m/t/1023182

There is a bug with DiscretizeGraphics where it isn't working on multi-polygons. Make this change in the definition of carto in either of the code blocks above and it should work

carto = DiscretizeGraphics[
   CountryData["France", {"Polygon", "Mercator"}] /. 
    Polygon[x : {{{_, _} ..} ..}] :> Polygon /@ x];